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Supersonic flow of a superfluid past a slender impenetrable macroscopic obstacle is studied in the 
framework of the two-dimensional defocusing nonlinear Schrodinger (NLS) equation. This problem 
is of fundamental importance as a dispersive analogue of the corresponding classical gas-dynamics 
problem. Assuming the oncoming flow speed sufficiently high, we asymptotically reduce the original 
boundary-value problem for a steady flow past a slender body to the one-dimensional dispersive 
piston problem described by the nonstationary NLS equation, in which the role of time is played by 
the stretched ^-coordinate and the piston motion curve is defined by the spatial body profile. Two 
steady oblique spatial dispersive shock waves (DSWs) spreading from the pointed ends of the body 
are generated in both half-planes. These are described analytically by constructing appropriate 
(-H . exact solutions of the Whitham modulation equations for the front DSW and by using a generalized 

Bohr-Sommerfeld quantization rule for the oblique dark soliton fan in the rear DSW. We propose 
an extension of the traditional modulation description of DSWs to include the linear "ship wave" 
£^ pattern forming outside the nonlinear modulation region of the front DSW. Our analytic results are 

supported by direct 2D unsteady numerical simulations and are relevant to recent experiments on 
Bose- Einstein condensates freely expanding past obstacles. 

C/3 ! PACS numbers: 47.40.Nm., 05.45.Yv, 47.40.Ki, 03.75.Kk, 03.75.Lm, 
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i. introduction 

In compressible fluid dynamics there are two canonical situations in which shock waves can be generated. In the 
| first case the formation of a shock occurs as a result of breaking of an evolving smooth or discontinuous profile of 
■ the density (velocity) and is described by the generalized solutions of the initial-value problems for the ideal fluid 
t^J- ' dynamics equations. The second type of shock waves occurs in a supersonic fluid flow past a body or as a result 
ON , of the motion of a piston within a tube filled with a liquid or a gas (see, e.g., [H, 0, and is associated with the 
boundary-value problems. In a viscous fluid, the shock wave can be represented as a narrow region within which strong 
dissipation processes take place and the thermodynamic and hydrodynamic parameters of the flow undergo a sharp 
\Q • change. However, if viscosity is negligibly small compared with dispersion effects, the shock singularity is resolved 
by a nonlinear wave train called a dispersive shock wave (DSW). A remarkable feature of the DSW is the generation 
ON . of solitons at one of its boundaries so that the whole structure can often be asymptotically described as a "soliton 
train" . 

An analytical theory of one-dimensional DSWs pioneered by Gurevich and Pitaevskii Q is based on the assumption 
that the oscillatory structure of a DSW can be asymptotically described by a modulated periodic (or, more generally, 
quasi-periodic) solution of the governing dispersive equation. The slow variations (modulations) of the travelling 
periodic wave parameters such as amplitude, wavenumber etc are governed by the so-called Whitham equations 
obtained by averaging of dispersive conservation laws over the period of the travelling wave. Analyzing the numerically 
observed structure of the dispersive shock wave, Gurevich and Pitaevskii proposed a special system of nonlinear free- 
boundary conditions for the KdV- Whitham system and obtained a global self-similar modulation solution for the 
problem of the decay of an initial discontinuity for the KdV equation. An analogous problem for the defocusing 
nonlinear Schrodinger (NLS) equation was formulated and solved in (see also a detailed analysis in 0, [f| where 

a different approach to the formulation of the step problem for the Whitham equations was used) . The modulation 
solutions describing more general cases of breaking of monotone and non-monotone initial profiles were obtained in 
StMELlKdV equation) and in [HCL1 (defocusing NLS equation) using Tsarev's generalized hodograph transform 
methodQ. 

The modulation theory of one-dimensional unsteady expanding DSWs proved to be very effective in different physical 
contexts ranging from shallow- water waves [l5[ to fibre optics to Bose-Einstein condensates (BECs) [H, E3 • in 
particular, it was successfully used in [l8| for the analytical description of the generation of dark solitons in quasi-lD 
transcritical BEC flows through wide penetrable potential barriers observed recently in the experiment [l||. 

The study of two-dimensional steady DSWs occurring in the supersonic dispersive flows past bodies was initiated 
in [20| where the stationary 2D system of the governing collisionless plasma equations was asymptotically reduced to 
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the ID KdV equation along the linear characteristics (Mach lines) with the stretched transverse coordinate playing 
the role of time (see also (21[) and then appropriate modulation solutions were constructed and interpreted in terms 
of the original steady 2D problem. 

Whilst an asymptotic description of supersonic dispersive flow past body in the framework of the weakly nonlinear 
KdV dynamics captures a number of essential features of the wave patterns arising in the flow, it may fail to provide 
an adequate description of the waves of sufficiently large amplitude. A different approximation not involving small- 
amplitude expansions, but instead, using the expansions in inverse Mach number, was proposed in |22j in the context 
of collisionless plasma dynamics. In 22] the problem of the supersonic dispersive flow past slender body was reduced 
to the so-called piston problem (this "hypersonic" transformation is known very well in classical gas dynamics — see, 
for instance, [2, 23]). In the Letter (24[ this transformation was applied to the problem of the supersonic 2D NLS flow 
past slender obstacle, which was translated into the piston problem for the ID dcfocusing NLS equation. It was also 
shown how the dispersive piston problem for the defocusing NLS equation can be asymptotically reduced to a much 
better understood initial-value problem. 

The present paper is devoted to a systematic study of the DSWs generated in the supersonic flow past extended 
bodies in the framework of two-dimensional defocusing NLS equation. The most relevant physical context of this 
problem is the description of the flows of Bose-Einstein condensates (BECs) past obstacles, which is currently a 
subject of intensive experimental and theoretical studies (see, for instance, fill [25], [26[ for recent experimental work 
and [3 H3, HH, H^, HO, HH H2, HH HH HH and references therein for some of the theoretical advances). It should 
be noted that the literature on this subject is growing too rapidly to reflect all recent advances. We also note that 
most of the existing theoretical work on the BEC flows past obstacles is concerned with the flows past small localized 
"impurities" with the dimensions of order of the healing length. In this paper, we consider an opposite situation, when 
the obstacle has the size much greater than the internal coherence length of the medium. This "slender body" problem 
is fundamentally important as a dispersive counterpart of the classical gas dynamics problem about the supersonic 
flow past a "wing" (see, for instance, [E0]) and has an advantage of the possibility of full analytical treatment. In 
addition, the solution of this problem elucidates the macroscopic mechanisms of the generation of dark solitons and 
"ship waves" in BECs observed in the numerical and physical experiments [25], [26], [2t| Hi], Hll HH HE] ■ Foreseeable 
direct physical applications could be connected with the BEC flows in atom-chip systems (see, e.g., (3(| H3 and 
references therein). 

The possibility of full analytical description of the 2D the supersonic NLS flow past body problem is based on the 
already mentioned "dispersive piston" approximation [24]. In the recent paper [381 ] . the dispersive piston problem 
for ID unsteady NLS flows was studied for the particular case of the piston moving with constant velocity (this 
corresponds to the flow past an infinite straight concave corner in the context of the present paper — see Section 
VI. A). In the present paper, full analytical modulation solutions will be constructed for this and other, more general, 
cases when the piston curve is a reasonably arbitrary unimodal function, which is necessary for the description of the 
supersonic NLS flow past a finite-length body. 

One of the unusual features of the NLS piston problem solution, not captured by the single- wave KdV approximation, 
is the generation of a nonmodulated nonlinear periodic wave in the region between the piston (body surface) and the 
trailing edge of the DSW for sufficiently large piston speeds. We show that this "transition wave" observed in the 
numerical solution in (38| actually occurs due to the reflection of a large-amplitude DSW from the piston surface — so 
that the interaction of the oncoming and reflected modulated waves necessarily leads to the formation of a region 
filled with purely periodic nonlinear oscillations. The occurrence of a nonmodulated nonlinear wave region in the 
similarity solutions of the defocusing NLS equation was first predicted in [6[ as one of the particular cases in the 
general classification of the decay of an initial discontinuity. 

In the recent studies [29, 31, 32, 33] of the supersonic BEC flow past localized obstacles two main distinct ingredients 
of the generated wave pattern have been identified and studied analytically and numerically: the so-called "ship waves" 
corresponding to the spatial Bogoliubov modes and generated outside the Mach cone, and oblique dark solitons 
generated inside the Mach cone and stretching behind the obstacle. An unexpected feature of these oblique dark 
solitons, established first numerically in [29[, is their apparent stability, in striking contrast with the well established 
notion of the absolute "snake" instability of two-dimensional dark NLS solitons leading to their decay into vortex- 
antivortex pairs [39], 0, Hl[ , • This apparent paradox was resolved in [3(| where it was shown that the presence 
of the background BEC flow with the velocity greater than certain "threshold" velocity stabilizes the dark soliton, so 
that it becomes only convectively unstable, i.e. practically stable in the reference frame attached to the obstacle. 

We note that in [2^, Hi], H2, HI] the ship waves and oblique dark solitons were studied as separate independent wave 
structures generated by an idealized obstacle of small size placed in the BEC flow. At the same time, the process of 
the generation of these wave structures, as well as the connection of their parameters with the geometry and size of 
the physical obstacle remained beyond the scope of the cited studies. In this paper, by considering an analytically 
tractable case of the supersonic NLS flow past a two-dimensional slender obstacle of finite size, we show that the ship 
waves and oblique dark solitons can be described as asymptotic far-field outcomes of the spatial "evolution" of two 
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separate DSWs spreading from the front and rear pointed ends of the body. In spite of their common origin, the front 
and rear DSWs evolve in drastically different ways: the front wave asymptotically transforms into a dispersing wave 
packet (effectively a "ship wave") while the rear one converts into a fan of dark solitons. This qualitative difference 
occurs owing to the fact that the front wave is developed from the compression "hump" forming due to the slowing 
of the oncoming flow near the increasing profile of the front part of the body while the rear wave evolves from the 
density dip forming behind the body. So in terms of the one-dimensional NLS equation, the front wave corresponds to 
the continuous spectrum of the associated Zakharov-Shabat linear spectral problem and the rear one — to the discrete 
spectrum. This is in striking contrast with classical dissipative gas dynamics where both shock waves spreading from 
the endpoints of the wing have essentially the same structure (see, for instance, [H, Hj]). 

We first develop the theory of the supersonic flow past a straight wedge by applying the similarity modulation 
solutions 0, 0] to the associated ID dispersive piston problem (see [H}). The comparison with full 2D numerical 
solution of the NLS equation with impenetrability boundary conditions at the body surface shows that the ID piston 
approximation describes the arising wave distribution remarkably well. 

Next we analyze the flow past a slender "wing" by constructing asymptotic ID analytical solutions for the front 
and rear DSWs and comparing them with the full 2D numerical solutions. We describe the front DSW behaviour by 
constructing an appropriate modulation solution of the dispersive piston problem with the piston curve corresponding 
to the body profile and then "translating" this solution in terms of the original 2D problem. 

The numerically observed wave distributions of the 2D NLS flow around the corner or the front edge of the wing, 
however, extend beyond the DSW region confined to certain boundaries, [y~ (x),y + (x)]. To describe the distribution 
of the wave crests outside the DSW region we extend the traditional Gurevich-Pitaevskii type formulation of the 
problem by complementing it by the modulation solution describing the distribution in the linear wave "packet" 
located outside the external DSW boundary y + (x). The lines of constant phase in this linear modulation solution 
determine the location of the small-amplitude wavecrests visible in numerical and physical experiments. Together 
with the DSW, they form a structure which eventually transforms into the universal Kelvin-Bogoliubov "ship wave" 
pattern [^3, [HJ • The far-field asymptotic behaviour of our nonlinear modulation solution describes the distributions 
of the wave amplitude in this "ship wave" as a function of the wing profile. 

Finally, we consider the rear DSW, which asymptotically decomposes into a fan of oblique dark solitons 29] . Instead 
of constructing the full modulation solution, we describe the asymptotic distribution of solitons in this fan using the 
generalized Bohr-Sommerfeld semi-classical quantization rule for the spectral eigenvalues obtained for the defocusing 
NLS equation in [131 |43j using the inverse scattering transform (1ST) formalism. 

Our analytical solutions are compared with the full numerical simulations of the 2D unsteady NLS flow past 
extended obstacles. 



II. FORMULATION OF THE PROBLEM 



We consider the supersonic NLS flow past an extended two-dimensional body with pointed ends (a "wing" ) . For 
simplicity we shall assume zero attack angle. 

We describe the flow dynamics by the multidimensional defocusing NLS equation in the canonical form 



iipt = -iA^ + MV 



(1) 



Since we shall be interested in the potential (vortex-free) flows it is convenient to transform Eq. |T|) to a hydrodynamic 
form by means of the substitutions 



tp(r, t) = y/n(r,t) exp (z9(r, t)) , 



u = ve. 



(2) 



where n(r,t) is the density of the "fluid" and u(r, t) denotes its potential velocity field, r = (x,y). We introduce 
normalized dependent variables h = n/riQ , u = u/c s , c s = ^/hq where no is the value of the density at infinity and c s 
is the corresponding sound speed. As a result, we obtain the system (we omit tildes for convenience of the notation) 



u t + (u ■ V)u + Vn + V 



m + V • (mi) = 
(Vn) 2 An' 
8n 2 4n 

V x u = 



0. 



(3) 



(here V = (d x ,d y )). 
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We assume the uniform oncoming flow with constant density n = 1 and the velocity u = (M, 0) directed parallel 
to x axis. Here M > 1 is the Mach number of the oncoming supersonic flow. The system © then should be solved 
with the boundary conditions at infinity, 



n — ► 1, u — > (M, 0) as 
and the impenetrability condition at the body surface S*: 

u-N| s = 0, 



(4) 



(5) 



where N denotes a unit vector of outer normal to the body surface. Similar to classical gas dynamics theory of 
supersonic flows (see 0] for instance) we shall be interested in an established, steady, wave pattern. Hence we confine 
ourselves to stationary solutions of the problem ©-([5]) and replace Eqs. ([3]) by their time-independent versions for 
n(x,y), u = (u(x,y),v(x,y)): 

(nu) x + (nv) y = 0, 



'.'7/7 



8n 2 



4n 



uv x + vv y + n y + 



8n 2 



= 0, 



0. 



(6) 



Uy - V x = . 



Let the shape of the body in the upper half-plane be given by a unimodal (one-hump) function: y = F(x) > for 
x S (0, L), F(0) = F{L) — and F(x) = for x ^ [0, L), L being the body length in dimensionless units (see Fig. 1). 
Thus we have N oc (F'(x), —1), and the boundary conditions ([!]), (O are transformed to 



u = M, 



at x 2 + y 2 



(7) 



v = uF'(x) at y = F(x). 
The flow in the lower half-plane can be considered independently in a completely analogous way. 



(8) 



u = (M,0) 



F(x) 



L x 



FIG. 1: Flow past a wing 



III. PISTON PROBLEM APPROXIMATION AND QUALITATIVE DESCRIPTION OF THE WAVE 

PATTERN 

The system ©-(H]) is still too complicated for a direct analytical treatment. However, when the flow can be con- 
sidered as highly supersonic the steady problem of the two-dimensional flow past slender body can be asymptotically 
transformed to a much simpler problem of ID "unsteady" flow along y axis with the scaled x coordinate playing the 
role of "time" [12]. To this end, we substitute into Eqs. ([5]) the new variables 



u = M + it! + 0(1/M), T = x/M, Y = y, 



(9) 
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assuming M Then to leading order we obtain 



VT + VVY +n Y+ [—-—) =0 



tit + (nv)y = 0, 
\ n YY \ _ n (10) 

Y 



Ul=0. (11) 

Equations (|10|) represent the hydrodynamic form of the ID defocusing NLS equation 

i^ T + i* yy - |^| 2 * = (12) 

for a complex field variable 



(13) 



and we can apply the well-developed analytical methods to its study. It is remarkable that in the case of a slender 
body, for which Ma = 0(1), where a = max|F'(a;)|, the boundary condition ([8]) reduces (to leading order in M _1 ) 
to the classical piston conditions (see for instance) 

v = Vp = df/dT at Y = f(T), (14) 

where the piston motion is described by the function /(T) = F(MT). Condition ((4]) at infinity transforms into 

n = l, v = as Y — > oo . (15) 

Thus, we have reduced the problem of the flow past slender body to the piston problem for ID flow along a tube 
with a piston moving inside it according to the law (I14| . In contrast to the classical gas dynamics, the piston problem 
is now posed for dispersive equations (|10p . 

The piston reduction for 2D hypersonic dispersive dissipationless flows was first introduced in [22| in a rather general 
form and in [24[ it was formulated in the present NLS context. In [38| the dispersive piston problem for ID defocusing 
NLS equation was studied for the simplest case of the piston moving with constant velocity (this corresponds to the 
flow past an infinite straight concave corner in the context of the present paper — see Section VI). In the subsequent 
sections, analytical modulation solutions will be constructed for this and other, more general, cases when the piston 
curve is a non-monotone function, which is necessary for the description of the supersonic NLS flow past a finite-length 
body (a "wing"). In classical viscous gas dynamics, the supersonic flow past a wing leads to the generation of two 




FIG. 2: Piston analogy in the problem of supersonic flow of dispersive fluid past body 

spatial shock waves (oblique jumps of compression) spreading from the front and the rear edges of the wing (see Q 
for instance) . In terms of the piston problem this corresponds to the formation of two shocks during two different 
phases of the piston motion: forward and reverse. 

Before we proceed with the quantitative analysis of this problem we briefly outline the qualitative structure of the 
dispersive flow past finite-length body using the theoretical results of [2(J HJ, [24[ ■ We assume that the length of the 
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body is much greater than typical dispersive (coherence) length of the medium. Then in dispersive hydrodynamics, 
both shocks spreading from the body edges resolve into expanding nonlinear oscillatory zones, the oblique spatial 
dispersive shock waves. At finite distances from the body surface these two spatial dispersive shocks have similar 
structure (see [H): each represents a modulated nonlinear wave acquiring a form close to a chain of oblique solitons at 
one edge of the oscillatory zone and degenerating into a linear wave at the opposite edge. However, as was indicated 
above, at large distances from the body the two dispersive shocks demonstrate drastically different behaviour: in the 
present case of the NLS hydrodynamics the dispersive shock spreading from the rear edge of the body transforms into 
the oblique soliton train while the dispersive shock forming at the front end of the body completely degenerates into 
a vanishing amplitude dispersing linear wave packet. 



IV. MODULATION THEORY FOR THE DEFOCUSING NLS EQUATION: ACCOUNT OF RESULTS 

The theory of DSWs is based on the study of a certain nonlinear free-boundary problem for the modulation 
(Whitham) equations — the so-called Gurevich-Pitaevskii problem. In this section we make a brief review of the 
relevant results of the modulation theory for the defocusing NLS equation which are necessary for the analysis of 
spatial DSWs generated in the steady supersonic NLS flow past slender body. A detailed derivation of the single- 
phase NLS modulation system can be found in [44| . 



A. Travelling wave solution and modulation equations 



The periodic travelling wave solution of the defocusing NLS equation (|10[) can be expressed in terms of the Jacobi 
elliptic function "sn" and is characterized by four constant parameters Ai < A2 < A3 < A4, 

n = j(A 4 - A 3 - A 2 + A1) 2 + (A 4 - A 3 )(A 2 - A x ) sn 2 (^J (A 4 — A 2 )(A 3 — Ai) 9, m) , (16) 



v = U--, (17) 

n 



where C = |(-Ai - A 2 + A 3 + A 4 )(-Ai + A 2 - A 3 + A 4 )(Ai - A 2 - A 3 + A 4 ), 



1 4 

= Y-UT-6 , U=-^\i, (18) 

i=i 



U being the phase velocity of the nonlinear wave and #0 initial phase. The modulus < m < 1 is defined as 

_ (A 2 - A0(A 4 - A 3 ) 



(19) 



(A 4 -A 2 )(A 3 -A 1 )' 
and the wave amplitude is 

a= (A4-A3XA2-A1). (20) 

The wavelength is equal to 

o_ 7 dX 2K(m) 

J ^(X-X 1 )(X-X 2 )(X-X 3 )(X 4 -X) v/(A 4 -A 2 )(A 3 - AO' 
A3 

K(m) being the complete elliptic integral of the first kind. 

In the limit asm->l (i.e. as A3 — > A 2 ) the travelling wave solution (I16p transforms into a dark soliton riding on a 
"pedestal" n : 

n = riQ ^ , (22) 

co$h 2 (^- s (Y-U s T-e )) 
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where the background density no, the soliton amplitude a s and the soliton velocity U s are expressed in terms of 
Ai, A 2 , A 4 as 



n 



r(A 4 -A 1 ) 2 



(A 4 - A 2 )(A 2 - Ai), U s = -(X 1 + 2X 2 + X i ). 



(23) 



Allowing the parameters Xj to be slowly varying functions of Y and T, one arrives at a modulated nonlinear periodic 
wave in which the evolution of Xj 's is governed by the Whitham modulation equations in the diagonal Ricmann form 



1,2,3,4, 



(24) 



which are obtained via the averaging of the NLS conservation laws over the period of the travelling wave solution 
([11]) (see [1, HH for the detailed description of the Whitham method) . The characteristic velocities can be calculated 
using the formula [9j, |4J| 



Vi(X) 

Substitution of Eq. into Eq. 



1 



di I U, i= 1,2,3,4, where di = d/8Xi. 



(25) 



gives the explicit expressions 

(A 4 -A 1 )(A 2 -A 1 )K 



Vi = 














*5> + 



(26) 



(A 4 -Ai)K-(A 4 -A 2 )E' 

(A 3 -A 2 )(A 2 -Ai)K 
(A 3 -A 2 )K-(A 3 -A 1 )E' 

(A 4 -A 3 )(A 3 -A 2 )K 
(A 3 -A 2 )K-(A 4 -A 2 )E' 

(A 4 -A 3 )(A4-Ai)K 
(A 4 -A 1 )K-(A 3 -A 1 )E' 

where E = E(m) is the complete elliptic integral of the second kind. The characteristic velocities (f26f are real for all 
values of the Riemann invariants therefore system (f2Tj) is hyperbolic. Moreover, it is not difficult to show using (|25|) 
that diVi > for all i so the NLS- Whitham system (|2"4")l . PB)) is genuinely nonlinear 

An asymptotic modulated wave solution is obtained by substituting the solution of the modulation equations (|24|) 
back into the travelling wave (fT6|) . We stress that initial phase 6*o in (fT8|) is "erased" in the averaging procedure so 
the resulting modulated wave is defined with the accuracy to an arbitrary shift within the wave spatial period. 

For the DSW analysis in the subsequent sections we shall need the reductions of the formulae (|26[) for the limiting 
cases when m = and m = 1. 

The harmonic limit m = can be achieved in one of two possible ways: one sets either A 2 = Ai or A 3 = A 4 . Then: 



When A 2 = Ai : V 2 = Vy = Ai + 



V a 



When A 3 = A 4 



V 3 = Vi - A 4 + 



Vi = 



A 3 + A 4 | 2(A 3 -A 1 )(A 4 -A 1 ) 
2 2Ai — A 3 — A 4 

3 1 3 1 

2^3+2^4, V" 4 = -A 4 + -A 3 

Ai + A 2 | 2(A 4 -A 2 )(A 4 -Ai) 
2 2A 4 — A 2 — Ai 



3 1 3 1 

-Ai + -A 2 , V 2 = -A 2 + -Ai 



(27) 



(28) 



In the soliton limit we have m = 1 . This can happen only if A 2 = A 3 , so we obtain: 

1 



When A 2 = A 3 : V 2 = V 3 
3 



-(Ai +2A 2 + A 4 ) 



(29) 



Vi 



1 



2 Ai + -A 4 , V 4 



3 x 1> 

2 A4+ 2 Al 
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Thus in both harmonic (m — > 0) and soliton (m — ► 1) limits the fourth-order modulation system (|24|) , (j26|) reduces 
to the system of three equations, two of which are decoupled. Moreover, one can see that in all considered limiting 
cases the decoupled equations agree with the dispersionless limit of the NLS equation (flQ|) . Indeed, the dispersionless 
limit of the NLS equation is the ideal shallow-water system 

Ut + (tiv)y =0, vt + vvy + fly = , (30) 
which can be represented in the diagonal form by introducing Riemann invariants 

A± = i»±Vn (31) 

§^ ± (A + ,A_)^0, (32) 

where 

V+ = ~A+ + \\- , V- = \\- + \\+ , (33) 



B. Hodograph transform and reduction to the Euler-Darboux-Poisson equation 

We fix two Riemann invariants, 

Ai = Aio = constant , A2 = A20 = constant (34) 
to reduce ([24]) to the system of two equations 

^ + K3(A 3 , A 4 )|i = , ^ + F 4 (A 3 , A 4 )|l = , (35) 

where Vs >4 (A 3 ,A 4 ) = V 3 ,4(Aio, A20, A3, A4). Applying the hodograph transform to system (|35|) one arrives at a linear 
system for F(A 3 ,A 4 ), T(A 3 ,A 4 ) 

dY dT dY dT 

A 3 ,A 4 )|- = 0, ^__y 3 (A 3 ,A 4 )^- = 0. (36) 
oA 3 0X3 0X4 0X4 

Now we make in (|36|) the change of variables 

Y-VjT = Wj, j = 3,4, (37) 
which reduces it to a symmetric system for Ws(Xs, A 4 ), W 4 (A 3 , A 4 ): 



Wi - W 3 Vi - V, 



; i,j = 3,4, Mi; ft = d/d\ . (38) 



The symmetry between Vj and Wj in (|38p and the "potential" structure (f2"5)) of the functions Vj implies the possibility 
of introducing a single scalar function <?(A 3 , A 4 ) instead of the vector (W 3 , W 4 ): 



or, which is the same, 



Wi=g + 2(Vi-U)^-. (40) 
OX, 



Then substituting (|25|) . (|40|) into (|38)) we arrive, taking into account (|21j) . at the Euler-Darboux-Poisson (EDP) 
equation for g(A 3 , A 4 ) first obtained in the present NLS context in Q (see also (T3 |) 
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The general solution of the EDP equation (|4Tj) can be represented in the form (see, for instance, [471 ]) 

" 3 <MA)rfA f 2 (A)riA 

V(A-A 3 )(A 4 -A) 7 V / (A-A 3 )(A 4 -A) ' 

where </>i,2(A) are arbitrary (generally, complex- valued) functions. 

As a matter of fact, the same construction can be realized for any pair of the Riemann invariants while the remaining 
two invariants are fixed. Moreover, equations ([37|) -(|41 |) turn out to be valid even when all four Riemann invariants 
vary 0, [l^]. This becomes possible for two reasons. Firstly, the NLS modulation system (fTTj) . (fl"4"|) is integrable 
via the generalized hodograph transform [l4[ which converts it into overdetermined consistent system (|38[) where 
i, j = 1,2,3,4, i ^ j. Secondly, the "potential" structure of the characteristic speeds (|25[) makes it possible to use the 
same substitution ([39]) for all i = 1,2,3,4 which results in the consistent system of six EDP equations (|41|) involving 
all pairs Aj, Xj, i ^ j. 

Thus, the problem of integration of the nonlinear Whitham system ([24)) with complicated coefficients (|26l) is 
essentially reduced to solving the classical linear EDP equation (|4"Tj) so practically one needs to express the functions 
0i,2 (A) in the general solution (j42|) in terms of the initial or boundary conditions for the NLS equation |T]). 

One should note, that classical hodograph solutions do not include the special family of the simple- wave solutions as 
the latter correspond to the vanishing of the Jacobian of the hodograph transform (Aj, Xj) i— ► (Y, T) (see, for instance, 
Q). However, the similarity solution can be formally included in the hodograph solutions in the generalized form 
(|37|) . Indeed, putting one of Wk — and setting constant all the Riemann invariants Xj with j ^ k one arrives at the 
similarity solution, in which Afe = Xk(Y/T) is implicitly specified by the equation Vu = Y/T. 



C. Free-boundary matching conditions for the modulation equations 



In the description of the DSW, the Whitham equations ([24)) must be equipped with certain boundary conditions 
for the Riemann invariants Xi These conditions are the NLS analogs of the Gurevich-Pitaevskii conditions [ij 
formulated for the KdV dispersive shock waves. To be specific, we formulate boundary conditions for the right- 
propagating DSW, which corresponds to the spatial DSW generated in the upper-half plane in the problem of the 
supersonic NLS flow past body. Without loss of generality we assume that the formation of the DSW starts at the 
origin of the (Y, T)-plane. In the Gurevich-Pitaevskii setting the upper (Y, T)-half plane is split into three regions 
(see Fig. 2): (-oo, Y~(T)), [Y~ (T),Y+(T)} and (Y+(T), +oo). 




FIG. 3: Splitting of the yT-plane in the Gurevich-Pitaevskii problem for the defocusing NLS equation. 

In the "outer" regions (—oo,Y~(T)) and (Y + (T),+oo) the flow is governed by the dispersionless limit of the 
NLS equation, i.e. by the shallow- water system ([32]) . ([33]) for the Riemann invariants A±. In the DSW region 
[Y~(T), Y + (T)) the averaged oscillatory flow is described by four Whitham equations (|24p for the Riemann invariants 
Aj with the following matching conditions at the trailing Y~{T) and leading Y + (T) edges of the DSW (see [1, fo r 
details): 



At Y = Y-(T): A 3 = A 2 , A 4 - A+, Ai = A_ , 
At y = F+(T): A 3 = A 4 , A 2 = A+, A x = A_ . 



(43) 
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Here X± (Y, T) are the Riemann invariants of the dispersionless limit of the NLS equation in the hydrodynamic form 
(f52| . ([55)) . The free boundaries F ± (T) are defined by the kinematic conditions 

dY~ dY + 

-^r = V 2 (Xi, A 2) A 2 , A 4 ) = V 3 {X U A 2 , A 2 , A 4 ) , = K 3 (A 1; A 2 , A 4 , A 4 ) = F 4 (A 1; A 2 , A 4 , A 4 ) (44) 

and so are the multiple characteristics of the Whitham system. The multiple characteristic velocities V2 = V3 and 
V3 = V4, in (j44f are explicitly given by ([29]) and (|28|) respectively. Determination of F ± (T) is an inherent part of the 
construction of the full modulation solution. We also emphasize that matching conditions (|4"3")) are consistent with the 
limiting structure of the Whitham system (j2~4")l at to = and m — 1 (see ([25)1 . ((23)) and reflect the spatial oscillatory 
structure of the DSW in the defocusing NLS hydrodynamics (as is known very well, such a DSW has a dark soliton 
(to = 1) at the trailing edge and degenerates into the vanishing amplitude harmonic wave (m = 0) at the leading 
edge— see 0,1, [H,p). 

One should mention that if one is interested only in the class of Y/T similarity modulation solutions arising in 
the decay of an initial discontinuity problem one can use, instead of (|43|) . a reformulation of the modulation problem 
as an initial-value problem for Xj, where three of the invariants are constant at T — and for the forth one the 
so-called "regularized" initial initial condition is used (see Hi Ell Ell ) ■ The resulting initial- value problem has the 
global expansion fan solution. This type of the problem formulation, however, seems to be less natural when one 
is interested in a more general (not self-similar) class of solutions when the integration of the modulation equations 
involves the hodograph transform (|37[) , (|38p (the poor compatibility of the initial- value problems with the hodograph 
method is known very well in classical hydrodynamics (see, for instance, [3j)). The free-boundary Gurevich-Pitaevskii 
type formulation (I43|) . on the contrary, is ideally compatible with the generalized hodograph transform as in any of 
the hodograph space coordinateplanes (A; , Xj ) it transforms into the classical Goursat-type characteristic boundary 
problem for the EDP equation [3, LL2J . 



V. ASYMPTOTIC REFORMULATION OF THE NLS PISTON PROBLEM AS AN INITIAL- VALUE 

PROBLEM 

The general dispersive piston problem (fl4|) , (|15p for the defocusing NLS equation (jT0|) is difficult to tackle directly. 
It is, therefore, desirable to reformulate it in terms of a much better explored initial-value problem. The key in this 
reformulation is the possibility to use the semi-classical Whitham description which is applicable when the character- 
istic piston displacements are much greater than unity while the piston speed is 0(1) (this formally corresponds to the 
supersonic flow past a slender body with the length L ^> M in our original setting formulated in Section II, however, 
one can expect that the results will be relevant to moderate body lengths as well). We now assume the qualitative 
picture of the flow described in the end of Section III and divide the upper part of the (Y, T)-plane in the piston 
problem into five distinct regions (see Fig. 4). In the regions 7 and V the flow is undisturbed so we have n = 1, v = 
there. The corresponding "dispersionless" Riemann invariants ([3Tj) are X± — ±1. In the region 777 for Y > f(T), 
the "gas" is put into "motion" by the "piston" moving according to Eq. (TlT| (we shall omit the the quotation marks 
for the terms related to unsteady gas flows henceforth) and near the piston the gas motion can be described by 
the dispersionless limit of the defocusing NLS equation (|32p . (f33|) . However, the formal solution of the nonlinear 
hydrodynamic-type equations (f3"2")l , (I33|) cannot be extended to the whole (Y, T)-plane because the ^-derivatives blow 
up along certain lines in this plane so the region 777, where the flow is smooth, is separated from the constant flow 
regions 7 and V by two DSW regions 77 and IV which spread from the points (0, 0) and (0, L/M), corresponding to 
the endpoints of the "wing" (strictly speaking, one should impose a certain restriction on the behaviour of f{T) near 
T = to have the front DSW emanating strictly from the point (0, 0) — this restriction will be explained in the end of 
this Subsection). The qualitative structure of these DSWs was described in Section II. We denote the leading (outer, 
i.e. facing the oncoming flow) and trailing (inner, i.e. facing the body surface) edges of the front DSW (region II) as 
Y^~(T) and Y7(T) respectively, and, similarly, for the rear DSW (region V) edges, we use the notations Y r ± (T). 

Now, the plan is to determine the flow parameters n p and v p at the piston surface and then to trace them back to 
T = using the solution of the dispersionless equations ([3"!?]) , (|3"3"|) . The kinematic condition (Ti"4"l) defines v p — df jdT 
so we just need to find the flow density at the piston. This can be done by considering the data transfer along 
the characteristics in the Gurevich-Pitaevskii setting of the problem where the entire wave pattern is asymptotically 
described by hyperbolic equations of hydrodynamic type (the NLS-Whitham system in the regions 77 and IV 
and the dispersionless limit of the NLS equation (f3"2")l in the regions 7, III and V). 
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We formulate the matching conditions for both DSWs using the general rule ((43)) . For the front DSW we have: 

At Y 



Similarly, for the rear DSW: 



Y f -(T) 
At Y = Y+(T) 

At Y = Y~(T) 



A3 — A2 , A4 — A+, Ai — A_ 
A3 = A 4 , A 2 = 1, Ai = — 1 . 

A3 = A2 , A4 — 1, Ai = — 1 , 



(45) 



(46) 



At Y 



Y+(T) 



A3 — A4 , A2 — A_|_, Ai 



A_ 



It then follows that to satisfy the governing equations (|24|) . (|32|) and the matching conditions (|45ll and (|46jl one has 
to put 



Ai = A. 



■1 



(47) 



within the respective domains of definitions of Ai (regions II and IV) and A_ (regions /, III and V). This condition 
(|47| of transfer of the Riemann invariant of the dispersionless system across the DSW replaces the traditional shock 
jump conditions for classical viscous shocks (see [49| for a detailed discussion of transition conditions across DSWs). 
Hence, we have at the "piston" v p /2 — ^JrTp~ = — 1 which yields the gas density 



17 CO, 



Y, + (T) 



L/M 



Yf(T) y;(t) 




FIG. 4: (Y,T)-plane of the NLS piston problem. Dashed line: the piston trajectory Y = f{T). The lines Yjr(T), Y±(T) are 
the edges of the front and rear dispersive shocks respectively. 



n p = (v p + 2) 2 /4 (48) 

in the region between the piston and DSW. Then using (|14p we get 

A_ = -1, A+ = df/dT +1 at Y = f(T). (49) 

We are now able to translate these boundary conditions at the "piston" into the equivalent initial conditions at T = 0. 
This problem for the system (|30|) can be easily solved using characteristics. Indeed, we have A_ = —1, hence A+ 
obeys the simple wave equation following from (J22J) (see, e.g., [IH) 



(50) 
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Solution of ([50)1 with boundary condition (|49|) is readily found using characteristics: 



A+ = /'(£) + 1 , Y= f(0 + + 1)(T - , 



(51) 



where £ is a parameter along the piston curve Y = f{T). Then, setting in (|51| T = we arrive at a parametric form 
of the equivalent initial distribution of the Riemann invariant A+ : 

A+(Y,0): A+ = /'(0 + l, Y = f(0-(lf(0 + l)Z. (52) 

Distribution (|52|) together with the initial condition 

A_(Y,0) = -1 (53) 

define, via (|3ip . initial conditions for the NLS equation in the hydrodynamic form (|10|) . It is important to emphasize 



Continuous spectrum: 
linear radiation 



Discrete spectrum 
dark solitons 




FIG. 5: Sketch of asymptotically equivalent initial conditions for A± in the problem of supersonic NLS flow past slender profile 

with /'(o) = .r (o = o. 

that initial conditions (|5^)) - ([55)) and the piston boundary conditions (fT4")l . (fTS"]) are equivalent only asymptotically, as 
our reformulation is made within the conditions of applicability of the Whitham modulation approach. One should 
also stress that in the context of the flow past body problem, the solution to the initial value problem (fTU|) . (|52 p - (|5"3"|) 
is defined only in the region Y > f(T), i.e. outside the body (piston). 

Now one can make some qualitative predictions about the asymptotic structure of the flow in the problem of the 
2D NLS flow past slender body. First we assume that /(0) = f(l) = /'(0) = f'(l) = 0, where I = L/M. Also to avoid 
unnecessary complications at this stage we assume that A+ (Y, 0) , specified parametrically by (|52")l , is a single- valued 
function (this restriction is not essential but it makes our analysis more transparent). Then it is not difficult to see 
from (|52p that the "translated" initial profile for the hydrodynamic Riemann invariant A + corresponding to the flow 
past a wing has the shape of a large-scale "bi-polar" pulse (see Fig. [SJ while the invariant A_ is constant. Note that 
the pulse is supported on the interval [—1, 1. T hen, the semi-classical approach to the inverse scattering transform for 
the defocusing NLS equation developed in |13l . [43| enables one to associate the "well" part of the initial profile of A + 
with certain distribution of dark solitons in the rear far-field asymptotic while the front "barrier" part is responsible 
for the linear dispersing radiation in the region II as T — > oo. The asymptotic formula for the amplitude distribution 
in the dark soliton fan generated out of the rear DSW will be presented in Section VII. The modulation solution for 
the front wave gradually transforming, via the nonlinear DSW stage, into the Kelvin-Bogoliubov "ship-wave" pattern 
(see HH, [3^) will be constructed in Section VI. 

For a more "realistic" wing shape (as in Fig. 1) we have /(0) = f(l) = f'(-0) = f'(l + 0) = but /'(+0) ^ 0, 
f'(l — 0) 7^ 0, where f'(a + 0) and f'(a — 0) denote the right and left derivatives of f(x) at x = a respectively. The 
qualitative behaviour of the solution remains the same but the the quantitative description undergoes some technical 
modification. Indeed, one can see from ({52")) that the discontinuity of the derivative /'(£) at £ = / implies that the 
rear endpoint of the wing maps back to an interval [Yi,} 7 !], where 



Y 2 = -l. 



Yi = -Gf'(! 



0) + l)l. 



(54) 



At these points the function A+ (V, 0) assumes the values 



A + (y 2 ,o) = i, A+(y 1 ,o) = i + /'(z-o) 



(55) 
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FIG. 6: Sketch of an asymptotically equivalent initial condition for A_| 
with /'(+0)/0,/'(Z- 0)^0. 



in the problem of supersonic NLS flow past slender wing 



On the interval [Yt, Y2] the function A+(Y, 0) is linear: 

A+(Y, 0) = 1 - | (1 + Y/l) for Y 2 <Y<Y U (56) 

and A+(Y2, 0) = 1 for Y < Y%. One can readily see that the function A+(Y, 0) is continuous everywhere except for the 
point Y — where the profile A+ (Y, 0) has a discontinuity: 

A+(0,0) = l + /'(+0) > 1, and X+(Y,0) = 1 for Y>0. (57) 

We also note that one can see from ([52")) that the point xq of the maximum of the body profile F(x) maps to the 
point Y"o = f(xo/M) — xq/M on the Y-axis so that \+(Yq, 0) = 1. A typical profile of the function A+(Y, 0) is shown 
in Fig. [5J 

For convenience of the presentation we shall assume that A+ ( Y, 0) is a single- valued function on the interval [Yi , 0] . 
This implies a restriction 

0< dA +( y '°) <oo for Y 1 <Y<0. (58) 
dY 

Then from ([52")) we obtain the corresponding condition for the function /(£) (i.e. for the body profile 

-AO 
i + + §£A0 



0< ; , j 777k , 3 ~„^ <°° fOT < ^ < Z . (59) 



One should stress, that actually there is no need for the function A+(Y, 0) to be one- valued as it is a formal projection, 
along the characteristics of the Riemann-Hopf equation ([50)1 , of the given physical distribution (Tl9"l) of A + specified on 
the piston curve. So our resulting formulae will not be restricted exclusively to the profiles satisfying the inequality 

We also formulate the condition necessary for the front DSW be generated exactly from the edge of the body at 
(0, 0). This is obtained from the condition that the profile A+(Y, T) breaks exactly at the initial moment T = as in 
Fig. 6. This is clearly the case if /'(+0) 7^ 0. However, if /'(+0) = then A+ can tend to unity at Y = according 
to the square root law, A + (Y) cx \[—Y as Y — > —0. This means that dY / d\ + \ T=a — at £ = and this condition 
can be satisfied if /"(£) — > 00 but /"(£)£ — > const where const can be equal to zero. In particular, such a behaviour 
takes place for /"(£) cx < j3 < 1. If dY/dA+| T=0 ^ then the wave breaking occurs at a later moment T = Tf, 
at the point Y = Y, which are determined by the equations 



dY 
dXl 



d 2 Y 

0, 



dA^ 



0. 



Simple calculation with the use of Eqs. ([ST]) yields the equation for £t>: 

(/'(&) + 2) /'"(&) = 4(/" (£ b )) 2 . 

When £b is found, then and Y> are calculated according to the formulae 

rp t 1 4/ /; (6) v _ ,^ , 2(3/(6) +2)/"fe) 

j6 -» + 3/'"(&)' 3/'«(&) 



In spatial x, y-terms this means means that the point of generation of the DSW is detached from the obstacle. 
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Y 1 = -l[l-\aM) , Y 2 = -l, ((,4) 



Example: parabolic profile — a "wing" 

We illustrate the described mapping of the body profile onto the initial profile of the Riemann invariant A+ by 
considering a parabolic "wing" with an opening angle a above the :r-axis and the length L so that the function F(x) 
in ([5]) is given by 

F(x) =ax(l-x/L), 0<x<L. (60) 

Then the piston function is f(T) = F(MT) = aMT{\ - T/l), where I = L/M. Now, the corresponding initial 
distributions of the "dispersionless" Riemann invariants A± are given by (|52")l , (f53"| . ([56]) . that is we have the following 
specification for A+ (Y, 0) 

A+(y,0): A+ = 1 + aM(l - 2£/t) , Y = -£ U + ^ (1 - 4£/[)\ for Y 1 < Y < , (61) 

A+(Y, 0) = 1 - ~ (1 + Y/l) for Y 2 < Y < Y l , (62) 
o 

X+(Y, 0) = 1 for - oo < Y < Y 2 and Y > . (63) 

Here (see ([55]). ([5i)» 

y; = -i (i - ^. 

so that the minimal value of A+ is A+(Yi, 0) = 1 — aM. Also we have 

r = / (2/2) - 1/2 = -2/2(1 - aM/2) , A+(r , 0) = 1 . (65) 

and A+(0, 0) = 1 + aM. We note that condition ([55]) is satisfied if the denominator in it is negative as £ — > 2 which 
implies a simple inequality aM < 2/7. 

VI. FLOW PAST STRAIGHT CORNER 
A. Analytical theory 

We first consider a model problem of the flow past an infinite straight corner specified by the function 

F(x) = 0, for x < 0; F(x) = ax for x > , (66) 

where a > is some constant (see Fig.[7])- To apply the piston approximation we need to assume that a ~ AT -1 <C 1 
so the piston curve in (|4"5|) is f(T) — aMT, i.e. the piston speed is 

v p = -± = aM, (67) 

Using (|48p we obtain that in the piston approximation the flow parameters in the region between the body surface 
and the DSW are simply 

U = M, v = v p = aM , n = n p = {Ma + 2) 2 /4. (68) 

Now, using Ip)2")) - (f53"]) we obtain the asymptotically equivalent initial conditions for the NLS equation (fT7j|) in terms 
of A± (see ([31])) 

T = 0: A_=-l; 

(69) 

A+ = A + > 1 for Y < , and A+ = 1, for Y > . 

Here 

A+ = 1 + aM . (70) 
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Dispersive shock 




FIG. 7: Supersonic dispersive flow past concave corner. The flow speed and density in the region between the corner and DSW 
are v p = aM, n p — (v p + 2) 2 /4. 



Of course, in the context of the flow past body problem, the solution is defined only for Y > Y p = aMT. Thus, the 
problem essentially reduces to the much studied problem of the decay of an initial discontinuity for the defocusing 
NLS equation (see 0,0]) with some restrictions for the domain of the solution. 
The relevant modulation solution has the form of a centered characteristic fan 

Ai=-1, A 2 = l, A 4 = A+, (71) 



or explicitly (see ([26|0 



where 



! = ^ 3 (-l,l,A 3 ,A+) (72) 



Y l a _ , An (l + aM-A 3 )(A 3 -l)K(m) 

T = 2 aM) (A3 - l)K(m) - aME(m) ' (?3) 



ro= 2 ^ i+ ° M -y, (74) 

aM(A 3 + l) V 7 



A + = A + < 3 4, 



" p < 2 
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FIG. 8: Behaviour of the Riemann invariants in the similarity DSW at some T > 0. The vertical double line at Y = Y p marks 
the "piston" (local body surface) position. Left: subcritical piston speed, v p < 2; Right: critical piston speed v p = 2 — formation 
of a "black" soliton at the trailing edge Y~ — Y p 

The DSW is confined to an expanding region t~T < Y < t + T, where the "speeds" (or rather slopes) r ± of the 
edges are calculated from (|75|) as the boundary values of the similarity variable r — Y/T: 



t = r(m = 1) = 1 + — , 



(75) 
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T ='■("» = °) = lTaM ' (76) 

We note that the trailing edge speed r~ is translated into the slope of the oblique dark soliton forming at the DSW 
edge facing the body surface 



= t~/M = 1/M + a/2. (77) 



(A 4 - A 3 )(A 2 - A x ) = 2(A+ - 1) = 2aM . (78) 



The amplitude of this soliton is 



The density profile in the soliton is defined by formula (|22p in which one substitutes the pedestal no = n p = 
(2 + aM) 2 /A, the amplitude a = 2aM and the 'velocity' U s — t~ . We note that it follows from (|77|) that in the 
flow past corner problem the oblique dark soliton is formed outside the "conventional" Mach cone defined by the 

50] where the 



slope l/VA/ 2 - 1 « 1/M. This is in an apparent contrast with the wave pattern described in [291 
oblique dark solitons were shown to be necessarily formed inside the Mach cone. However, in 23 . 50l] the oblique 



dark solitons were considered to be generated by the point-like obstacle. In that case the background flow density 
and, correspondingly, the sound speed were equal to unity so that the Mach number in the background flow was 
everywhere equal M . In the present case of the flow past corner, the oblique dark soliton is generated on a non-unity 
background which results in a different value of the local sound speed and, therefore, in the changed definition of the 
Mach cone which is now specified by the local Mach number Mi = Mj Jn^,. As a result, the adjusted Mach angle 

becomes 1/ \/M? — 1 w 1/M + a/2 which coincides with the oblique soliton slope 1(77)1 . Thus, in the supersonic NLS 
flow past a corner an oblique dark soliton is formed along the actual Mach line. This agrees with the result in Q 
where it was shown that the trailing dark soliton in the DSW moves with the sound speed. Since s + = t + /M > s~ 
the implication of this fact is that the DSW is located entirely outside the Mach cone. 

The schematic behavior of the Riemann invariants in the modulation solution (|7ip . ()T3)1 is shown in Fig. 8. As was 
already mentioned, in the context of the supersonic flow past body (or the piston) problem the solution (|73|) is defined 
only for Y > v p T — aMT. Then from the condition r _ = aM we obtain the critical value v p = aM = 2 for which 
the greatest dark soliton in the DSW is generated right at the body surface (see Fig. 8, right panel). Incidentally, 
this value of v p also implies that the density at the minimum of this greatest soliton turns zero which constitutes the 
appearance of a vacuum point at the trailing edge of the DSW 6]. Indeed, from (|16|) . the minimum of the density 
in the travelling wave solution is n m i n = ^(A 4 — A3 — A2 + A1) 2 . Substituting (|71[) we obtain the distribution for the 
local minima of n in the DSW, 

n min (m) = i(aM - A 3 (m) - l) 2 , (79) 

where the dependence A3 (to) is given by (|T4[) . Then the requirement that n m i n (l) — immediately yields aM = 2. 
Generally, setting in ((79)) aM > 2 one gets from n m i n (m) = 

< 1 , (80) 



(aMf 

i.e. the vacuum point occurs inside the DSW. Since at the vacuum point we have Ai = — 1, A 2 = 1, A 3 = aM — 1, 
A4 = 1 + aM the phase velocity ((18)) at the vacuum point is U* = aM, i.e. is equal to the piston velocity. This seems 
to imply that the DSW gets attached to the piston and is realized only partially with the modulus ranging from to 
in*. However, it turns out that one cannot attach the partial DSW directly to the piston, instead, one should introduce 
an additional periodic "transition wave" with m = to* between the DSW and the piston (38|. As a matter of fact, the 
vacuum point is present at each period of this transition wave. The generation of a nonmodulated periodic nonlinear 
wave in the piston problem can be explained in the following way. The vacuum point phase velocity U* coincides with 
the nonlinear group velocity only when to* = 1, i.e. for the dark soliton, when the multiple characteristic velocity 
V2 = V3 = U is nothing but the soliton speed. Generally, for to* ^ 1 one has V2 7^ V3 7^ v p and, therefore, one should 
introduce a reflected wave. As a result of the DSW reflection from the piston (body surface) one generally would 
get a two-phase wave region characterized by six Riemann invariants (see, for instance [g, |48| | for the corresponding 
Whitham equations), with two of them changing. However, the requirement of self-similarity of the modulation 
solution in the problem of the supersonic flow past straight corner imposes the restriction that only one Riemann 
invariant can change. This implies that the two varying Riemann invariants in the general two-phase modulated 
solution must coincide with each other with the consequence that there is only one oscillating phase described by 
four distinct constant Riemann invariants (the varying multiple Riemann invariant can be ignored as it essentially 
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FIG. 9: Riemann invariants for the supercritical "piston velocity", v p > 2. The periodic transition wave occupies the region 
[Y P ,Y*] where Y* = r*T (see (EE])). 



describes the propagation of the vanishing amplitude linear wave packet against the cnoidal wave background). So, 
as a result of nonlinear wave interaction one effectively gets a non-modulated finite-amplitude periodic wave, which 
in the present context can be viewed as a nonlinear standing wave. The behaviour of the Riemann invariants in 
the described "supercritical" modulation solution is schematically shown in Fig. 9. The region of the intermediate 
periodic wave expands with the speed 

'• = «-',U»-M + cM) = .M - (aM !^ K ^ > E( ,„. ) , (81) 

where to* = 4/(aM) 2 (see (l80l) ). The transition wave amplitude is (see (1201) ) 

a* = (A 4 -A 3 )(A 2 -A 1 )=4, (82) 

and it does not depend on the value of aM > 2. The latter only affects the transition wave width r*T (and the local 
wave shape, via to*). This is in striking contrast with the classical dissipative piston problem, where the density jump 
across the shock increases without limitations as the piston velocity grows. We also note that in the physical xy-plane 
the transition wave is located between the body surface and the centered line with the slope s* = t*/M. An explicit 
expression for the oscillating density profile in the transition wave in physical xy plane is 

n = 4sn 2 (aM(y - ax), m*) , (83) 

and the wavelength of the transition wave (83)) in any x-section is calculated as £ (see (|2"Tj) ) and is given by 

£ = ^M K{m * ] ■ (84) 

In conclusion we note that the described transition wave solution actually represents part of the special similarity 
solution obtained in as a particular case in the decay of an initial discontinuity problem (see the case 6 in the full 
classification of [6]). The transition wave was also very recently observed in [38] in the numerical simulations of the 
dispersive piston problem. 



B. Comparison with numerical solutions 



We have performed two series of numerical simulations. Firstly, we constructed the full unsteady numerical solutions 
for the 2D NLS flow past corner for M — 10 and different values of the corner angle a. Secondly, we performed parallel 
numerical simulations of the associated ID piston problem (fT2"j) - (|15p for the corresponding values of the piston velocity 
v p = Ma. In both cases we have used finite-difference codes with the impenetrability condition ijj = at the body 
(piston) surface. The results of the numerical simulations have been then compared with the analytical modulated 
solution obtained in the previous subsection. We note that from the numerical point of view it is more convenient to 
perform simulations for a symmetric wedge with the opening angle a above the rr-axis and to use the wave pattern in 
the upper half-plane for the comparison. 
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FIG. 10: Dependence of the flow density n p near the corner surface on the vertical velocity component v p for the flow with 
M — 10. Dashed line: dependence (|68|) obtained in the dispersive piston approximation. Circles: data obtained from the full 
2D numerical solution. Solid line: n p (v p ) - dependence for the classical dissipative piston problem 



First, we have made a comparison for the DSW transition condition, which in our case is expressed by the formula 
(|68|) specifying the parameters of the constant flow between the DSW and the corner surface provided the oncoming 
flow parameters are u = AI, v = 0, n = 1. The comparison of the dependence n p (v p ) for M = 10 with the numerical 
data for the density in the 2D flow past a wedge is shown in Fig. 1101 This is also compared with the dependence n p (v p ) 
for the classical dissipative piston problem for polytropic gas with 7 = 2, corresponding to the dispersionless limit of 

the NLS equation. The corresponding classical piston jump condition is specified by the equation v p — (n p — 1)^/-^-^- 

The comparison is made for M = 10. One can see excellent agreement between the analytical dispersive piston curve 
and the numerical data obtained from the full 2D simulations of the flow past corner problem. At the same time one 
can see noticeable departure of the dispersive piston curve from the classical piston curve. The numerical simulations 
data and the dispersive piston curve split at v p = 2 (i.e. at a = v p /M = 0.2, which also agrees with our solution as 
for Ma > 2 the theory predicts the formation of a transition wave so that the region of a constant flow between the 
corner and DSW disappears. 
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FIG. 11: Left (color online): 2D density plot for of the supersonic (M = 10) NLS flow past a wedge with the opening angle 
above the a;-axis a = 0.1. Right: density profile n(y) at x — 50. Top: analytical modulated solution; middle: numerical solution 
of the associated ID piston problem; bottom: ID cut of the full 2D solution at t = 15. Points y~ and y + mark the boundaries 
of the DSW predicted by the modulation solution. The body surface (piston) is located at about y — 5. 



In Figs. [TTrfT3l the 2D density plots (left) and ID cross-section density profiles (right) are presented for the flows 
with M = 10 past corners with a = 0.1,0.2 and 0.3 respectively. The analytical solutions (top panel) are compared 
with the numerical solution of the asymptotically equivalent ID dispersive piston problem (middle panel) and with the 
x-section of full 2D solution. One can see that ID numerical dispersive piston solutions agree remarkably well with the 
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FIG. 12: Left (color online): 2D density plot for the supersonic (M = 10) NLS flow past a wedge with the opening angle above 
the :r-axis a — 0.2. Right: density profile n(y) at x = 30; top: analytical modulated solution; middle: numerical solution of 
the associated ID piston problem; bottom: ID cut of the full 2D solution at t — 15. Points y~ and y + mark the boundaries of 
the DSW predicted by the modulation solution. The body surface (piston) is located at about y = 6. 




FIG. 13: Left (color online): 2D density plot for the supersonic M — 10 NLS flow past a wedge with the opening angle above 
the x-axis a — 0.3. The dashed line marks the end of the transition wave predicted by the theory. Right: density profile n(y) 
at x = 30; top: analytical modulated solution; middle: numerical solution of the associated ID piston problem; bottom: ID 
cut of the full 2D solution at t — 15. Points y~ and y + mark the boundaries of the DSW specified by the modulation solution. 
The body surface (piston) is located at about y = 9. 

results of the full 2D simulations. The agreement between the analytical solutions and numerical simulations is also 
very good in the DSW region (we note that the exact position of the wave in the analytical solution is determined up 
to a characteristic coherence length (soliton half- width) as the initial phase 6>o in (| 16[) is not denned by the modulation 
theory). The predicted occurrence of the vacuum point at the body surface at aM = 2 and the generation of the 
non-modulated transition wave for aM > 2 are seen very well in Figs. [T^l EH The predicted position y* = t*x/M of 
the right boundary of the transition wave (see (|8ip) also agrees very well with the numerical simulations — see Fig. 1131 

Next, in Fig. [14] the comparisons for the amplitude (|78p and slope (|77|) of the first dark soliton in the DSW as 
functions of the corner angle are presented. One can see that the agreement is excellent for the soliton amplitude 
and quite good for the slope. One should emphasize that the accuracy inherent the hypersonic approximation (flU|) 
implies that the amplitude formula (|75|) is defined with the accuracy 0(1/M) while for the slope s~ given by ([77]) 
the accuracy is 0(1/M 2 ). Since the slope formula is s~ = 1/M + a/2 one can expect that a noticeable discrepancy 
between analytical and numerical values for of may be the case for small angles, say for a < 0.1. Of course, this 
will contribute, on level of 0(x* /M 2 ), to the error in the analytical determination the spatial y-location of the oblique 
dark soliton at some x-cross-section made at x = x* (as in Figs. [TTHT3")) . We note that the analytically predicted 
soliton location is also subject to an arbitrary, up to a typical wavelength, shift inherent in the modulation theory. 

One should also note some important feature of the wave pattern that is not captured by the modulated solutions as 



20 




a a 

FIG. 14: Parameters of the first dark soliton in the DSW as functions of the corner angle. Left: the soliton amplitude a~; 
Right: the soliton slope s~ . The numerical values are taken at x — 50 

seen in the right upper panels of Figs. [TTI fLll Indeed, one can see noticeable small-amplitude oscillations beyond the 
outer harmonic edge y + of the DSW (as defined by the modulation theory). In the theory of one-dimensional DSWs 
these linear oscillations are usually ignored. However, in the considered here 2D problem these linear oscillations 
represent an essential part of the observable wave pattern (see the left panels in Figs. [TTr[T3)) and should be taken into 
account. A similar wave distribution was considered recently in [3l|, [32|, [33[ in connection with the Bogoliubov-Kelvin 
"ship waves" generated by a point-like obstacle placed in the supersonic BEC flow (see also in [26| the discussion 
of the experimentally observed patterns). An extended modulation solution describing the combined wave pattern 
including both the DSW and the linear "ship-wave" distribution will be constructed in the next section. 

It is worth noting that in the strongly nonlinear region near the the wedge boundary at large x one can see the 
oscillations of the dark soliton crest lines (see the density plot in Fig. [T2] (left panel)). This is the manifestation of 
the so-called "snake" instability of dark solitons with respect to bending disturbances [H, 0, 5l[ . However, for large 
enough oncoming flow velocity these unstable disturbances are convected by the flow along solitons and hence they 
become just convectively unstable in the reference frame related with the obstacle [3(3]. Therefore, for the considered 
here large Mach numbers, the DSW structure can be regarded as effectively stable and thus can be treated as a 
modulated stationary solution of the 2D NLS equation. 

VII. FLOW PAST WING 




FIG. 15: (color online) Supersonic NLS flow past a wing: density plot. The oncoming (from the left) flow speed is M = 10. 
The green dashed line shows the wing profile specified by Eq. (160[) 

Now we consider supersonic flow past an extended slender finite body — a wing. From the very beginning we assume 
zero attack angle so without loss of generality we shall consider the wave pattern only in the upper half-plane. The 



21 



density plot for supersonic (M = 10) flow past the wing having a symmetric parabolic form specified by the function 
(JHOJ) with a — 0.15 and L = 100 is shown in Fig. [T5] One can see that the wave pattern agrees with the qualitative 
predictions made in Section V using the inverse scattering transform reasoning applied to the asymptotically equivalent 
initial data of the type shown in Fig. [5l Indeed, one can see the front DSW, similar to that in the straight wedge 
case described in Section VI and the fan of oblique dark solitons spreading from the rear edge of the wing. Unlike 
the straight wedge case, though, the front DSW is not characterized by a constant jump of density n and velocity 
v across it, so the depth of the oscillations decreases as the distance from the generation point at (0.0) increases. 
As a result, the front wave degenerates into a small-amplitude disper sing wave with the distribution of wavecrest 
wave having the form similar to the "ship- wave" pattern described in [3ll |32| . The length of the wing used in the 
simulations is not sufficiently large to identify the details of the intermediate front and rear DSWs. However, we shall 
construct full modulation solution for the front DSW and, by considering its asymptotic behavior for large x,y will 
derive the amplitude and wavelength distributions applicable to the ship-wave pattern. For the rear DSW, instead of 
constructing full modulation solution, we shall take advantage of the semi-classical Bohr-Sommerfeld type distribution 
fl3l l43l ] for the distribution of eigenvalues in the Zakharov-Shabat scattering problem. 

The crucial difference between our consideration in this paper and the results obtained in earlier papers [Sana ib 
on dark solitons and ship waves is that here we asymptotically solve the boundary value problem for the 2D NLS 
equation and express the parameters of the resulting wave distributions in terms of the initial profile, while the 
previous papers were concerned with the study of certain particular solutions of the 2D NLS equation. 

A. Flow past front edge of a wing 

1. Formulation of the problem 

To model the flow of a superfluid past the front edge of a wing we consider the function F(x) of the type shown 
in Fig. 16 (left) so that F = for x < 0; F'(+0) = a < 1, F'(x) > for < x < x and F'(x) = for x > x . 
Then one can readily see from (|52[) that for highly supersonic flows the asymptotically equivalent initial condition 
for A + = |d + tJu has the shape shown in Fig. 16 (right). The other Riemann invariant A_ = —1 (see (|53|) ). Also 
Yo = /(Co) ^ Co where Co = xq/M. In terms of the piston problem this corresponds to the forward motion of the 
piston. The initial piston velocity is v p = Ma (as in the problem of the flow past straight corner with the angle a) 
but then the motion of the piston slows down until it eventually stops at T — Co ■ To avoid unnecessary complications 
connected with the formation of the transition wave we shall assume that aM < 2. 




Y ¥ 

FIG. 16: Left: front edge of a wing in an upper half-plane. Right: asymptotically equivalent initial condition for A+. 

As was explained in Section V, it is clear from the IST-based reasoning that the disturbance caused by the front 
edge of the wing in the supersonic NLS flow will eventually (for T 1) transform into a linear dispersive wave 
radiation. However, for an intermediate values of T the spatial "evolution" of this disturbance leads to the formation 
of a DSW having the structure similar to that generated in the flow past straight corner described in the previous 
Section. Thus, remarkably, even in this "solitonless" configuration, the DSW and dark solitons still form, albeit as an 
intermediate wave pattern. While in the evolutionary problems this wave pattern is transient, in our 2D stationary 
problem the intermediate front DSW exists for all times and transforms into linear waves only at large distances from 
the body. The essential difference is that, due to the presence of the spatial scale xq, the corresponding modulation 
dynamics is no longer self-similar resulting in the wave parameter variations along the wave crest lines which now 
have curved geometry. In particular, one can expect that the oblique dark soliton forming at the trailing edge of the 
DSW will initially (i.e. at x = 0,y = 0) have in the Whitham approximation the slope s~ = a/2 + 1/M and the 
amplitude a~ = 2aM as in the corresponding straight corner case, but as the distance from the body increases, its 
amplitude and slope will both decrease and asymptotically one can expect that a~ — > and s~ — > X/M as x — > 00. 
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2. Modulation solution 



We use (|45l) . (|47|) to formulate matching conditions for the Riemann invariants in the front DSW (we shall omit 
the subscript in Yj): 

At Y = Y-(T): A 3 = A 2 , A 4 = A+, A x = -1, , 
At Y = Y+(T): A 3 = A 4 , A 2 = 1, Ai = -1, lSOj 

where A+(Y, T) is the solution of the simple- wave equation (150)) with the initial condition A+(Y, 0)) defined by (I52p . 
Thus we have for A + an implicit representation 

Y-i(3A+-l)T = ™(A+), (86) 

where w(X + ) is the inverse function to \+(Y, 0) (note that for general non-monotone initial profile A+(Y, 0) one would 
need to consider two monotone branches of w(\+) separately, but in our case of the profile shown in Fig. 16, right, 
there is only one branch specified by (|52p ). Schematic behavior of the Riemann invariants corresponding to the 
matching conditions (|85[) is shown in Fig. 17. 



A 




FIG. 17: Schematic behavior of Riemann invariants in the modulation solution for the front DSW. 

Importantly, the modulation problem ()24l) , (|97| is no longer self-similar so one should use the hodograph transform 
to solve it (see Section IV. B). First, since Ai = — 1 and A 2 = 1 satisfy both modulation equations (|24|) and matching 
conditions (|85[) we have Ai = — 1 and A 2 = 1 everywhere, hence there are only two modulation equations for A3 and 
A4 left to solve. These transform via the substitution (see (13710 

Y - V 3 (— 1, 1, A 3 , A 4 )T — JY 3 (A 3 , A 4 ) , Y- Y 4 (-l, 1, A 3 , A 4 )T = 1Y 4 (A 3 , A 4 ) (87) 

into a system of two linear partial differential equations for 1Y 3i4 (A 3 , A 4 ) 

1 dW 3 _ 1 dV 3 1 cW 4 _ 1 dVj 

Wi - W 3 d\i Vi - V 3 8X4 ' W 3 -Wi dX 3 V 3 -V 4 dX 3 ' [ ' 

The boundary conditions for equations (|88p are obtained by considering the hodograph solution (|87[) at the free 
boundaries Y^ and applying to it the matching conditions At the trailing edge Y = Y~ (T) we have A 3 = 1 and 
Y 4 (— 1, 1, 1, A 4 ) = |(3A 4 — 1) (see (|2T)|) ') so that the second equation ([57)) becomes 

Y — i(3A 4 - l)T = W 4 (l, A 4 ) . (89) 

Comparing (|59")1 with the simple-wave solution ([5d| at Y = Y _ , where A 4 = A + (see ([55]) ) we obtain the boundary 
condition for (|88|) 

Wi{l,Xi)=w{X A ). (90) 
At the leading edge Y = Y+ we have A 3 = A 4 and (see (l28l) ) 

Y 3 (-l, 1, A 4 , A 4 ) = Y 4 (-l, 1, A 4 , A 4 ) = 2A 4 - 1/A 4 = Y*(A 4 ) . (91) 
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The multiple characteristic velocity V*(\4) determines the speed of the leading edge (see (04])). Note that since 
A4 > 1 we always have d±V* > and for the initial data of the type shown in Fig. 16 (right) the speed of the 
leading edge at T — is V*(A + ) which is the greatest characteristic speed in the system. Therefore the characteristic 
dY/dT = V*(A + ) is not intersected by other characteristics of the family dY/dT — V4 so the equation of the leading 
(harmonic, m = 0) edge of the DSW is simply 

Y- (2A+ - 1/A+)T = 0. (92) 

Substituting A + = 1 + aM we obtain the slope of the outer (facing the oncoming flow) edge of the DSW in the 
physical x, y plane as 

+ _ y _ 2(aM) 2 +4(aM) + l 
S x ~ M(l + aM) ■ [ ' 

Thus, the outer edge of the spatial DSW is determined by the opening angle a alone and does not depend on the 
specific body contour (indeed, (|93[) coincides with the expression for the slope of the outer edge of the DSW generated 
by the flow past straight infinite corner with the angle a (see (|76jl ). 

The obtained leading edge equation ([92]) should be consistent with the hodograph solution ([87]) considered at m = 0. 
Then comparing (JH7J) for A 3 = A 4 = A + with (J52J) we get 

W 3 {A+,A+) = W 4 (A+,A+) = 0. (94) 



Equations (|90|) and (|94|) provide boundary conditions for the linear system (|88|). Using transformation (|40|) . which in 
our case is explicitly represented as 

^(A 3 ,A4) = < ? + 2[y l (-l,l,A 3 ,A4)-i(A 3 + A 4 )]^- ) i = 3,4, (95) 

the system ({88} is further reduced to the EDP equation ([4T|) for the potential function g(A 3 ,A4) (see Section IV. B 
for details). Now we need to translate boundary conditions (|90|) and (|94|) for the hodograph equations (|88|) into the 
boundary conditions for the EDP equation. 
Substituting ([95]) into ([90]) we obtain 

g(l, A 4 ) + 2(A 4 - 1) = ^(A 4 ). (96) 

Ordinary differential equation (|96p is readily integrated to give 

where we have chosen the constant of integration such that g(l,A + ) = (this requirement is not essential). Now, 
without loss of generality we take the general solution (|42]) of the EDP equation ([41]) in an equivalent form 

5(A 3 , A 4 ) 



0i(A)dA A f <^ 2 (A)dA 



V(A 4 -A)(A 3 -A) J N /(A-A 4 )(A-A 3 ) ' 
± A4 

where ^1,2 (A) are arbitrary (generally, complex- valued) functions. From (|94|) we obtain </>i(A) = 0. Next, applying 
boundary condition ([97]) we arrive at the integral Abel equation (see, for instance, [5l|) for <j>2 (A) 



A+ A+ 

<t> 2 {\)d\ _ 1 f w{z) 



y/(X- A 4 )(A- 1) 2y/\4^1 J VI^T 



dz . (99) 



Ai 



The solution to equation (]99p (obtained via the inverse Abel transform for fa (A)/vA— T) is 



1 /" — wfz) . 
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Substituting <fii = and 4>2(X) given by (jlOOp into the general solution (|98|) and changing the order of integration we 
obtain a compact representation for the solution to the EDP equation for the problem of the flow past front part of 
the wing 

g(Aa,M) = — 7= / , . K — — —t I dz , (101) 

ttVA 4 - 1 J Vz-X 3 \{z - A 3 )(A4 - I) J 

where K(z) is the complete elliptic integral of the first kind. Now, formulae (|101[) . ([55| and ([57)) provide the exact 
implicit modulation solution to the NLS initial value problem with the initial profile of the type shown in Fig. 16 
(right). Strictly speaking, one should now show that the obtained solution is global, i.e. the mapping (A3, A4) 1— > (Y, T) 
specified by ([57]), and (|10ip is invertible for all T. However, instead of giving full mathematical proof of the 
invertibility of the hodograph transform (|87[) for our solution, it seems to be more instructive just to show that 
the obtained modulation solution has a physically meaningful asymptotic behaviour for T » 1, which, apart from 
providing us with the useful information about distributions of physical parameters at large distances from the body, 
will be a convincing enough indication that the solution is valid for all T . 

To study the long-time behaviour of the obtained solution we express T from the hodograph formulae (|87[) . ([55]) as 

T _ Ws - Wa JV 3 - |(Ag + A 4 )]^ - [V, 1(A 3 + A 4 )]^ 

where we have denoted Vj = Vj(—1, 1, A3, A4), Wj = Wj{Xj,, A4) for brevity. Next, substituting the solution (|10ip into 
(|102p we obtain an explicit expression for T in terms of A3, A4. Analysis of this expression shows that T — > 00 implies 
A3 — > A 4 . Since the wave amplitude a = 2(A4 — A3) and the modulus m = 2a /[(X4 — 1)(A3 + 1)] (see ftT§\> ) we obtain 
that a ^ 0, m — > as T - > 00 everywhere except for a small vicinity of the trailing edge point where A3 — > 1, so one 
has a — > but m — * 1. That means that the front DSW asymptotically transforms into a vanishing amplitude linear 
wave packet (the asymptotic behaviour of the trailing soliton will be considered separately). 

Indeed, a straightforward analysis shows that for the obtained solution 1^3,4^3, Xi)/T — * as A 3 — > A4 (i.e. 
T — » 00). Then we have from the hodograph solution (|87p to leading order in 1/T 

T»l: ^ = ^ 3 (-l,l,A 4 ,A 4 )=2A4--^. (103) 
J A4 

Next, expanding (|102p for small A4 — A3 -C 1 we obtain, after some algebra, the leading order asymptotic behaviour 
(provided A4 is not too close to 1) 

a=^MXi), (104) 

where 




Asymptotic behaviour (|103p . (|104|) is consistent with the modulation theory for linear waves (see for instance [3J]). 
Indeed, using the definitions of the phase velocity U (fTT)) and the wavenumber k = 2ir/£ where £ is the wavelength 
pTjl one can see that in the linear limit A3 — > A4 one has 



2tt 



A 3 = A 4 : fc = — — — =2^1-1, so f/ = A 4 = v/l + P/4, (106) 

i^^— 1, 1, A4, A4J v 



the latter being the linear dispersion relation of the NLS equation (fT2|) . u>o(k) = kU = fcyl + /c 2 /4. Then the right- 
hand side of (fT03]) . 2A 4 - 1/A 4 = (1 + k 2 /2)/^/l + h 2 /A, is nothing but the linear group velocity u' (k) so (fT03l) is 
simply the similarity solution of the kinematic modulation equation kx + u)' (k)kY = for the linear wave packet. 

The asymptotic behaviour (|104|) of the amplitude is also consistent with the linear wave energy conservation 
law dra 2 + 9y(wp(fc)a 2 ) = 0. However, the function A(X4) defining the relation of the asymptotic wave amplitude 
distribution with the body profile cannot be determined within the linear theory and requires the full nonlinear analysis 
presented here. One should mention that the eventual transformation of the front DSW into a linear radiation 
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also agrees with the general reasoning of the inverse scattering transform method as the initial conditions of the 
type described in the beginning of this section (see also Fig. 10) represent a "solitonless potential" having only 
a continuous spectral component. Indeed, formulae ()103|1(|105[> could also be obtained via the inverse scattering 
transform formalism but the employed here method via the solution of the Whitham equations appears to be more 
direct and efficient for the purpose. 

Finally, using (I103j) - (|106|) we represent the asymptotic amplitude and wavenumber distributions explicitly in terms 
of the original spatial variables x, y to perform later a comparison with the numerical simulations of the 2D NLS flow 
past slender obstacle: 



x,y > 1 : 




k 



8 -16, 



where r = M— . (107) 

x 



One can see that k — > as r — > 1, the latter being the Mach line in the hypersonic approximation. This will also 
emerge in the next section where the trailing (soliton) edge of the front DSW will be shown to asymptotically approach 
the Mach line as x — > oo. A remarkable feature of the asymptotic wavenumber k distribution in (|107p is that it does 
not depend on the shape and size of the body (provided the conditions of applicability of the piston approximation 
are satisfied). This will allow us to construct an analytic description of the universal "ship- wave" pattern generated 
in the supersonic NLS flow past slender bodies. 





FIG. 18: Comparisons for the asymptotic ((x/M) S> 1) amplitude a (left) and wavelength 2n/k (right) distributions in the 
front DSW generated by the wing with a parabolic profile (|60[) (L = 100, a = 0.15) placed in the supersonic NLS flow with the 
Mach number M — 10. The comparisons between the asymptotic modulation solution (|107[) (solid line) and numerical solution 
(circles) are made for a fixed x = 50. 

The amplitude distribution a(x,y) in ()107|) . on the contrary, depends, via the function A(Xi), on the wing profile. 
We stress that in spite of the fact that the asymptotic distribution a(x,y) satisfies the amplitude equation of the 
linear modulation theory (see [H), the determination of the amplitude dependence on the boundary conditions (i.e. 
the determination of the function A(Xi)) has required full nonlinear modulation analysis. To explicitly evaluate the 
function A{\^) for the parabolic profile (|60p we just need to know the function w{z) entering the integral in (|105[) . 
Since w(z) is the inverse of X+(Y, 0) on the interval [lo,0] it is readily obtained as 

W(Z) = ~2^M {aM + 1 - z )( z -^f)- ( 108 ) 

We recall that I — L/M, where L is the length of the "wing". Then the integral in (|105|) is evaluated explicitly to 
give: 



^' = V^^ < A «- 1 »' ,4 (Ti( 1 + ° M - A «) 3/2 < 8A '- M, + 2 >) ■ (109) 

The comparisons of the asymptotic distributions (|107j) , (|109j) with the distributions of a and k obtained from the 2D 
numerical solution are shown in Fig. 1181 One can see a very good agreement for both distributions. 



3. Trailing edge 

The leading (outer) edge y + (x) of the DSW is determined by formula To complete the modulation solution 

we need to determine the trailing (inner) edge y~(x) defined by the soliton condition m = 1. As in the straight corner 
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case, we shall mainly be concerned with the amplitude of the trailing soliton and its slope as its actual position might 
differ significantly from the curve y~ (x) obtained from the modulation theory due to the loss of the initial phase (see 
the relevant discussion and comparisons with the numerical solution in Section V.B). 

In the piston problem terms, we are going to find the curve Y = Y~(T), where A3 = A2 = 1. Remarkably, the 
equation for Y = Y~(T) and, as a result, the parameters of the trailing soliton can be found directly, without using 
the full modulation solution obtained in the previous subsection. We use the fact that in the front DSW one has 
A2 = 1 (see the previous subsection) so the matching condition (|85l) at the trailing edge assumes the form: 

At y = Y~(T) : A 3 = A 2 = 1, A 4 = A+, Ai = -1, (110) 

where \+(Y,T) obeys the simple- wave equation (|86p . On the other hand, the curve Y = Y~(T) is specified by the 
kinematic condition (f44|) in which we set the values of A/s from (|110|> . As a result we get a closed system 

rlY 1 

Y-±(3\ + -l)T = w(\ + ), _ = -(i + A + ) (111) 

along the trailing edge. We introduce in Qlll|) Y = Y S (X*), T = T S (A*), A + = A*, where A* is the parameter along 
the trailing edge curve so that Y S (A + ) — 0, T S (A + ) — (since A + (0,0) = A + — see Fig. 14). Next, eliminating Y' a we 
obtain a single ordinary differential equation for T S (A*) 

(A*-i)r; + |r s + w '(A*) = o, T s (A+) = 0, (112) 
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which is readily integrated to give 



A+ 



Ts = j£ 1)3/2 J(z-l)^w'{z)dz. (113) 

A* 

Next, substituting ()113|) into the first equation (|111|1 we obtain the function io(A*) in the form 

Y s = ~(3\*-l)T s (\*)+w(\*). (114) 

Thus, equations (JTT3J)— (JTHJ) specify the DSW trailing edge {Y = Y~(T): Y = Y S {X*),T = T S (A*)}. Correspondingly, 
the geometric location y~{x) of this edge in the physical x,y-plane is given by 

y = y-(x): y = Y s (\*), x = MT s (\*). (115) 

Within the NLS modulation theory the position of the trailing edge determines, up to an inherent phase shift, the 
location of the trailing dark soliton. Thus, equations (|113p - (|115p define the geometric shape of this spatial trailing 
dark soliton. We note that, unlike recently found oblique dark solitons generated in the 2D supersonic NLS flow past 
small obstacles [29j], [3(|, and stretching along straight lines, the trailing dark soliton in the front DSW has a curved 
contour in the xy-plane. In fact, the "bending" of this 2D soliton has the same nature as the speed variations of 
a ID soliton propagating through a non-uniform medium. Here the non-uniformity is due to the large-scale density 
variations in the flow past extended obstacle. 

Using (|113p we obtain an implicit expression for the variations of the trailing dark soliton amplitude a~ = 2 (A* — 1) 
along the wave crest line y~ (x) 



1+aM 

2 3 / 2 M 



x = 



(a-) 3 / 2 

1+0-/2 



{z-lf/ 2 w'{z)dz. (116) 



The relationship between the local slope s of the trailing dark soliton and its amplitude is given by (see (|llip ) 

•--^-bC +•-/«>■ < 117 > 

Since close to the the origin, (x, y) — > (0,0), we have A3 — > A 2 = 1, A 4 — > A + = 1 + aM we get for the soliton 
amplitude a~(0, 0) = 2(A 4 — A3) = 2aM, so the initial slope of the trailing edge is s~(0) = 1/M + a/2, i.e. it coincides 
with the slope of the dark soliton in the DSW generated in the flow past straight corner (see ((77)) as one can expect 
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(note that this result does not have much practical significance as the modulation theory performs rather poorly for 
small x, y). 

Next, since the integral in (|116p is 0(1) we conclude that a~ ~ x~ 2 ^ 3 — > as x — > oo i.e. the trailing dark 
soliton amplitude vanishes along the line y~ (x) while its slope asymptotically approaches the Mach line of the highly 
supersonic undisturbed flow: y~ — > x/M as x — > 1. 

There still remains an issue of the transition from the amplitude decay a ~ x~ x l 2 (|107p for the major part of 
the the DSW to the decay a ~ x~ 2 / 3 (|116p for the trailing dark soliton at trailing edge. This matching requires a 
detailed analysis of the asymptotic behaviour of the hodograph solution (|87p in the small vicinity of the singular point 
A 3 = A4 = 1. Such an analysis, whilst being relatively straightforward, is beyond the scope of the present paper. 

For the parabolic wing profile (|60p the function y~(x) defining the location of the trailing dark soliton is given by 
(|115p . where for T S (X*) we obtain from (|113[) by using formula l|108p for the inverse function w{z): 



T, 



I 



30oM 



(10 - 3aM) 



aM 
X* - 1 



3/2 



12A* + 1MM + 2 



(118) 



and for Y S (X*) we have (TlMf . 

And, finally, for the trailing soliton amplitude a~(x) at j; > 1 we obtain from Q116P (or directly from (jTTSJ)) a 
simple implicit formula 



I (, ,{2aM\ 3/2 

(10-3aA/)(_) -6a- 



IbaM - 10 



(119) 



One can see that at x = one has a = 2aM and a — > as x — > 00 as predicted by the general theory. In particular 
for x 3> 1 we have the asymptotic behaviour of the amplitude along the soliton wavecrest, 

Z(10-3aAf) V /3 2aM 

30a J xV* 1 UJ 

The comparison of the amplitude behaviour (|119p along the DSW trailing edge with the numerical simulation data 




10 20 30 40 50 60 70 



FIG. 19: Comparison for the amplitude decay along the trailing dark soliton. Solid line: asymptotic modulation solution (|119p . 
circles: the amplitude values from the direct numerical simulation. 

is shown in Fig. 1191 One can see a good agreement between the analytical curve and numerical solutions for large 
enough (as expected from the range of validity of Eq. (|119p ) values of x. 

We note that dependencies (|116p and (| 1 1 7[) have been derived under an implicit assumption that the "time" x/M 
of the establishment of the trailing dark soliton in the DSW is much less than the typical modulation "time" scale 
~ L/M, i.e. under the assumption that the DSW is fully established. This assumption works quite well for the 
straight wedge type profiles studied in Section VI but it may fail for the wing-type profiles with sufficiently rapidly 
decaying derivative so that condition 1 is n °t satisfied for a significant part of the profile (or, in terms of 

asymptotically equivalent initial conditions (|52|) the inequality being \d\+(Y,0)/dY\ < 1). In that case, the trailing 
soliton establishes itself very slowly and realizes only asymptotically for x 3> L (see [52j for the analysis of a similar 
issue in the context of the KdV equation). Taking into account that the amplitude of this "slowly developing soliton" 
decreases with x on the scale ~ L, its behaviour for finite x/L could actually be rather well approximated by the 
linear theory (see the next section). At the same time, one should stress that the dependence (| 1 1 6[) of the values of 
the trailing soliton amplitude on the obstacle size and shape cannot be found within the linear theory. Determination 
of this dependence requires full nonlinear analysis (either modulation or IST-based) - see the discussion in section 
VII A.2. 
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4- Extension of the modulation solution: the "ship wave" pattern. 

The modulation solution obtained in section VII A. 2 is denned within the domain y~(x) <y< y + (x) and implies 
that the wave amplitude vanishes at the outer (leading) DSW edge y + (x) and outside of the DSW region flow is 
assumed to be constant. At the same time, the boundary y = y + (x), associated with the linear group velocity, is not 
a wavecrest line so it is clear that one should be able to extend the wave crests beyond the DSW boundary. Indeed, 
it is clearly seen from the results of numerical simulations (see the density plots in Figs. 11-13) that the wavecrests 
do not stop at the external boundary of the modulation solution and the small (but quite noticeable) oscillations are 
present outside the DSW. To resolve this apparent contradiction one can notice that the vanishing of the amplitude 
at y — y + (x) for the DSW modulation solution does not necessarily imply that the actual wave amplitude turns into 
zero, this simply means that the oscillations are linear. To capture these linear oscillations occurring for y > y + (x) 
we introduce a small-amplitude wave packet as a natural extension of the DSW, and will use the linear modulation 
theory for its description. 

In linear modulation theory the equation for the wave amplitude is decoupled from the equation for the wavenumber 
(see [H) so one can put a = and consider the "wave conservation" law separately. We note that such an extension, 
whilst being automatically consistent with the DSW modulation solution at y = y + (x), is not quite trivial as the 
linear modulation theory is not valid inside the DSW region, even in a small neighborhood of the zero-amplitude 
leading edge y + (x) — see [53J. We note that the modulation solution for the wavenumber in the linear wave packet 
has already been obtained, this is equation (jl03|) (see the explanation after formula (105)), so we simply postulate 
that this solution describes the wave distribution for y > y + (x). 

In effect, modulation solution (|103[) enables one to derive the two-dimensional "ship-wave" pattern generated by 
the front edge of the obstacle. To this end, we notice that, up to an arbitrary initial phase Go G [0, 2ir], the local 
angular phase of the two-dimensional "travelling" wave is given by (see (|T8"|) ). 



x 



(-) = M = /•■„ ( // - Vjg) =kyy-^l + %-x, (121) 




hence the wave vector of the modulated linear wave is equal to 

• "" ' ' (122) 

As in the 2D theory of ship waves produced by a point-like obstacle in the supersonic NLS flow [3l|, [32| , we introduce 
the angle x between the radius- vector r and the x axis, i.e., the flow direction, and the angle n between the wave- vector 
k and ~x axis (see Fig. 20, left)): 

r = (r cosx, r sinx), k = (— |k| cost?, |k| sin 77). (123) 
Then Eq. (|122j) leads to the following expression for the wave vector length 



|k| = 2 ^ /M2 , COt2 ^ 1 (124) 
sin 77 

in the hypersonic approximation. One should emphasize that the wavenumber k defined by Eq. (|106p and occuring 
in the zero-amplitude limit (|103[) of the one-dimensional piston approximation of the DSW modulation solution is 
consistent with the y-componcnt k y of the full two-dimensional vector k (| 1 22[) . Hence the substitution of 



A 4 = ^1 + fc2/4 = yjl + |k| 2 sin 2 77/4 = McotTj (125) 

into Eq. (102) yields the relationship between x an d 77 

tanx = 2 cot 77 - tan^/M 2 . (126) 



Now we notice that Eqs. (|124[) and (|126[) are nothing but the highly supersonic approximation of the "ship wave" 
theory developed earlier [3l|, [32j for the case of a localized point-like obstacle. Indeed, in this theory the length of the 
wave vector is given by the expression 



cos 2 r/- 1 (127) 
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FIG. 20: Left: The wave crest geometry in the NLS ship wave pattern. The wave vector k is normal to the wave crest line 
which is shown schematically by a curve. Right: the theoretical wave crest lines in the ship- wave pattern (solid lines); the DSW 
boundaries determined by the modulation solution are shown by the dashed line. 



which in our hypersonic approximation can be easily transformed to 

|k| = 2sin7 ? y / M 2 cot 2 rj - 1. (128) 
This expression is approximately equal to (|124[) if sin r/ = 1 , that is 

tan?7>l, or \k y /k x \ > 1. (129) 

In a similar way, the relation 

tony- fl + M a /2)fa»"* M30) 

taUX - Af 2 -(l + |k| 2 /2) { ' 

between the angles x an d f] for the point-like obstacle case in the hypersonic limit can be cast into the form 



2 tan r\ tan 3 r\ 

tanv = k — = -, (131) 

A tan 2 77-I M 2 (tan 2 ? ? - 1)' v ' 



and again this formula is reduced to (|126[) under condition f| 129[) which means that the flow parameters change 
much slower in the x-direction than in the y-direction what is assumed in our approach. Thus, we have arrived at 
a remarkable result: the solution of the Whitham equations describing the DSW region, turns out to coincide, for 
i,!/> 1, with the corresponding approximation of the linear ship wave theory describing the waves outside the DSW. 
Thus, the far-field asymptotic solution (| 1 0T|) is not restricted to the DSW region and can be used for the description 
of the whole flow at the distances sufficiently far from the front edge of the wing. 

This observation permits us to extend the the wave crest lines to the whole region outside the Mach cone. It remains 
only to show that the ship wave pattern produced by a slender body can be approximated by the pattern produced 
by a point-like obstacle. To this end we turn to the formula for oscillations of density in the ship-wave theory (see 
Eq. (20) in H|), 

f ]/(k)|k|V k "- dk 

071 J (k-U) 2 -|k| 2 (l + |k| 2 /4)(2^) 2 1 ) 

where U = (M, 0) and 

V(k) = / V(r)e- lk r dr (133) 



is the Fourier image of the potential V(r) created by the obstacle. Integral over wave vector length |k| can be estimated 
as contribution of the poles in (|132p which depends on the dispersion relation only. Moreover, for obstacles with a 
sharp form their Fourier images must include wide range of harmonics and hence they are smooth functions of k. 
Therefore in the integration over directions of k performed for |k ■ r| 3> 1 by the stationary phase method (see [32| ) 
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the main contribution is given by a stationary point of the phase k • r which, again, does not depend on the function 
V(k). This yields the relation (|130[) and, subsequently, the parametric formulae for the wave crest lines, 



V 



46 

w 

46 



cos 77(1 — M 2 cos 2-q), 
sim](2M 2 cos 2 77 — 1), 



where 6 = 2tt, 47r, 67r, . . ., the wavevector |k| is given by Eq. (|127[) and 7/ changes formally in the range 

— arccos(l/M) < f] < arccos(l/M). 
In the limit cos 77 3> 1 /M the crest lines take a parabolic form 



x{y) 



6 

2M 



M 
26 



y 



(134) 



(135) 



(136) 



This limit corresponds to the region located not too far from the front edge of the obstacle. 

In the opposite limit, when cos 77 — ► 1/M the crest lines converge asymptotically to the straight lines parallel to the 
Mach cone lines: 



x 

M' 



(137) 




FIG. 21: (color online) Comparison of the universal "ship-wave" pattern (|134[) for a slender body (dashed lines) with the 
wavecrest lines in the supersonic NLS flow past the front edge of the parabolic wing profile (|60|l with a = 0.15, L = 100. The 
oncoming flow speed is M — 10. 



It is important to note that the obtained expressions for the geometry of the wave crest lines do not depend on 
the opening angle a of the obstacle and, hence, on the slope of the outer edge y + (x) of the DSW. Indeed, as we 
have shown, although in the Whitham approximation the amplitude of the wave vanishes at y + (a;), the curves of the 
wave crest lines can be continued outside this DSW boundary where they represent the spatial distribution of the 
small-amplitude linear waves. For this reason the linear ship wave pattern can be viewed as a natural continuation of 
the wave crest distribution in the oblique spatial DSW. However, when we go inside the DSW region along the wave 
crest line, the amplitude of the DSW gradually increases and the linear approximation loses its applicability. In an 
established (strongly nonlinear) DSW, the shape of the wave crests is determined by the shape of the obstacle (see 
section VI. A. 3). But for the profiles with sufficiently rapidly decaying derivative the DSW establishment "time" x/L 
is rather large (see the discussion in the end of section VI. A. 3) so the linear approximation works quite well in a wide 
region around the front tip of the obstacle including the neighborhood of the outer boundary of the DSW. This is 
illustrated in Fig. [5T] by the comparison of the analytical predictions given by Eq. (|134|) with the results of full 2D 
numerical simulations. 
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B. Flow past rear edge of a wing 



Now we consider the DSW generated by the flow past the rear edge of the wing. The corresponding initial profile 
of the Riemann invariant A + is given by Eq. ([52"]) for Y\ < Y < Yq and by (JHU) for Y% < Y < Y\. We recall that 
Y = f(x ) - x , Yx = -(§/'(/ - 0) + 1)1, Y 2 = -I and A+(Yi,0) = 1 + f'(l - 0) (see ([55])). The second invariant 
is constant, A_ = — 1. A typical form of the function X+(Y, 0) is schematically shown in Fig. [17] (right). The evolution 
of the "potential well" A + (Y, T) leads to the wave breaking at 



Y±-Y 2 



K+Cl.-lJ-V^A+Cn.O),-!) 



= /, Y b = Y 2 + V+{\,-l)T b = {). 



(138) 



which simply means that the rear DSW spreads directly from the rear endpoint of the wing (see Fig. [4J . A typical 
profile of A+ (Y, T b ) is shown in Fig. [22] (left). 
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FIG. 22: Left: profile of the rear edge of a wing in the upper half-plane. Right: asymptotically equivalent initial condition for 
A+. 



Thus, for T > T b one has a DSW forming behind the body (see Region IV in Fig.fJJ. This DSW can be described 
by the modulated travelling wave solution analogous to that constructed in the previous section for the front DSW. 
The main difference is that now one has the Riemann invariants A3 and A2 varying within the modulation solution 
while Xi = — 1, A 4 = 1 (see Fig. [T5] (right)). As a result, the modulation solution yields that as T — ► 00 one has 
A2 — * A3 (i.e. m — ► 1) everywhere except some small vicinity of the leading edge where A3 = 1 and m — > 0. That 
means that the rear DSW for 1 » 1, j » 1 asymptotically transforms into a soliton train (a fan of oblique dark 
solitons). Of course, such a behaviour is to be expected as the initial profile of A+ (see Fig. |2"21 fright)) corresponds to 
a large-scale 'potential well' in the associated scattering problem in the Zakharov-Shabat 1ST formalism for the ID 
NLS equation, and, therefore, leads to a semi-classical distribution of the bound states, each linked to a dark soliton 
in the NLS equation solution [l3|, [43| . 




FIG. 23: Left: profile of the Riemann invariants X± at the point of wave breaking, T = Tb. Right: schematic behaviour of the 
Riemann invariants in the modulation solution for rear DSW, T > Tb 



Thus, if one is interested in the asymptotic structure of the flow in the region far enough from the body where the 
rear DSW transforms into a "fan" of spatial solitons well separated from each other, there is no need to derive the full 
modulation solution. As was shown in Refs. [l^ill^l, each soliton in the soliton train evolving from the initial "well" 
is parameterized by the eigenvalue A = A^ found from the generalized Bohr-Sommerfeld quantization rule, consistent 
with the Whitham approximation used before, 



V^A - A+)(A - A_) dY = 2ir(k + ±), k = 0,l,...,K, 



(139) 
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where in our case A + = X+(Y, 0) is given by (|52]) , (I56)) , A^ = — 1, and the integration is taken over the cycle around 
two turning points defined by A = A + (Y, 0) . The fc-th soliton amplitude a k is related with the eigenvalue Afc by 

o fc = 1 - X 2 k (140) 

Returning to the spatial coordinates we find the profile of the Afc-soliton in the train as (see [l3j]) 

1 - A 2 

n k (x, y) = 1 5 — fc , (141) 

C osh 2 [JT=>? k {y-{X k /M)x)y 

that is the "fan" of spatial dark solitons is made of soliton "feathers" lying asymptotically along the lines 

y = (\ k /M)x, k = 0,1, . . . ,K, (142) 

in the upper half-plane and symmetric "fan" of solitons is generated in the lower half-plane. 

Remarkably, the distribution (|139[) is invariant with respect to the evolution, up to a breaking point at T = Tb, of 
the Riemann invariant A + , described by the simple- wave equation (|50[) (which is consistent with the dispersionless 
limit of the NLS equation (jT2J)) . Indeed, it is not difficult to show that (|50|) implies that 

^/v / (A-A + )(A + l)dr = 0. (143) 

The property (|143|) can be viewed as a semi-classical analog of isospectrality of the ID NLS evolution (see [43l|). 
Thus, the initial profile of A+ is defined up to the deformation ([50]) and, thus, should not necessarily be a single- 
valued function as in Fig. [53] (see also the discussion in Section V). 

For the parabolic profile ([60]) . the function X+(Y, 0) corresponding to the rear part of the wing is specified by 
formulae (|61|). (|62|) in the interval Y% < Y < Yq and A + = 1 outside of this interval. It has its minimum at Y\\ 
A+(Yi,0) = 1-aM. 

Now the integral in (|139[) is evaluated in a closed form giving the equation for the bound states A = Afc.: 



— 4?(A-l + aM) 3/2 (3aM + 8A + 2) = 27r(fc + i), fc = 0, 1, . . . , N. (144) 
loaM 

The physically meaningful roots Afc lie in the interval 1 — aM < X k < 1 . It immediately follows from the requirement 
A+ > A_ = — 1 that one must also impose a restriction that aM < 2 (for aM > 2 the description should be modified 
as the vacuum point appears at y = 0). The greatest root A at has the value close to unity so that the number of 
solitons in the fan can be estimated by putting Afc = 1, fc = N in (1144|) . i.e 

N « — — (aM) 1 / 2 (10 + 3aM) (145) 
2n 15 

The semi-classical formula (|144|) . strictly speaking, is asymptotically valid as long as N ^> 1, which, by (|145p presumes 
rough general criterion l(aM) 1 / 2 3> 1. However, as is often the case with the Bohr-Sommerfeld type distributions, 
formula (I144[) works reasonably well for a much broader range of parameters. Say, for I = 10, a = 0.15, M = 10 one 
has just 3 physical roots of the equation (|144p . which agrees with three dark solitons observed in numerical solution 
(see Fig. IT51) . The comparisons between the predictions of p44[) and the numerical simulations data for the amplitudes 
a k and slopes s k of the oblique dark solitons are presented in the Table below. 



fc 


A fc 


dfe = 1 - 




a k (num) 


s/s = Afc /A/ 


Sfc (num) 





0.2915 


0.9150 




0.9170 


0.0291 


0.02 


1 


0.7101 


0.4957 




0.5689 


0.0710 


0.06 


2 


0.9649 


0.0688 




0.1903 


0.0964 


0.09 



One can see that, taking into account the inherent in the hypersonic approximation error 0(1/ M) for the soliton 
amplitude and 0(1/M 2 ) for the slope, the comparison should be viewed as quite favorable. 
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VIII. DISCUSSION 

In this paper, we have constructed an asymptotic theory of the supersonic flow of a supcrfluid past slender bodies. 
The theory is constructed in the framework of the 2D defocusing NLS equation with the impenetrability condition at 
the body surface and the condition of an equilibrium steady flow with Mach number M at infinity. The description is 
made under the following assumptions: M>1, a<l, Ma = O(l), where a is the body slenderness parameter (e.g. 
the opening angle of a "wing" or a wedge) . Under these assumptions we have asymptotically (with respect to the small 
parameter 1/M) reduced the original two-dimensional stationary boundary- value problem for the time- independent 
2D NLS equation in the x, j/-plane with the oncoming flow along the x-axis to the dispersive piston problem for 
ID defocusing NLS equation, in which the role of time is played by the stretched ^-coordinate, T — x/M, and the 
spatial variable is the transverse coordinate y. The flow is globally described using the semi-classical approximation 
of the NLS equation, when the solution is governed by the dispersionless limit equations (the shallow-water system) 
in the regions of non-oscillating flow and by the Whitham modulation equations in the regions of dispersive shock 
waves, representing rapidly oscillating expanding nonlinear wave structures. We use the so-called Gurevich-Pitaevskii 
formulation of the problem to match the solutions of the Whitham equations with the solutions of the shallow-water 
equations at free boundaries. The full modulation solutions are constructed and analyzed for two canonical cases 
of the supersonic flow past bodies: the flow past infinite straight corner (a wedge) and the flow past a wing. Our 
analytical solutions are supported by direct 2D unsteady numerical simulations. 

We now summarize the main results of the paper. 

• We have shown that the highly supersonic NLS flow past 2D slender bodies is accompanied by the generation 
of two DSWs with contrasting asymptotic properties. 

• By making the comparisons of the numerical solutions for the 2D problem of supersonic flow past infinite wedge 
with the ID numerical and analytical modulation solutions of the associated dispersive piston problem we have 
shown that the piston problem describes the arising 2D wave patterns remarkably well for sufficiently large 
Mach numbers. 

• Using the dispersive piston approximation, we have constructed exact modulation solutions for the problems of 
the supersonic NLS flow past a straight infinite wedge and a slender "wing" . 

• By analyzing the asymptotic behaviour of the obtained modulation solution for the front DSW in the flow past 
a wing we have derived the distributions of the amplitude a and the wavenumber k far enough from the front 
edge of the wing (Eq. (|W|l ): 

x,y>l: a = (^pj ' A^ T + V J 2+ -^j , fc^I^( T + Vt 2 + 8) 2 -16, where r = M~ , (146) 

where the function A(£) is given by Eq. (|107p . These distributions describe the Kelvin-Bogoliubov "ship wave" 
pattern and relate it, via the function A(£), with the geometric parameters of the wing 

• The distribution of oblique dark solitons in the rear DSW is obtained using the generalized semi-classical Bohr- 
Sommerfeld quantization rule. 

The developed in this paper theory could find the applications to the description of Bose-Einstein condensates 
behaviour in current experiments on loading of ultracold quantum gases in traps, their coherent manipulation and 
transport. Such processes are now under intense investigations in atom chips — microfabricated, integrated devices 
in which electric, magnetic and optical fields can confine, control and manipulate cold atoms. Understanding of the 
interplay of dispersive and nonlinear properties in Bose-Einstein condensate dynamics is of crucial importance for the 
effective use of these devices which have very promising technological applications. 
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